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Abstract 

Computational  design  of  new  materials  relies  on  accurate  descriptions  of  interatomic 
potentials.  Such  potentials  can  be  realized  within  the  Multi-Body  Expansion  (MBE) 
framework,  where  the  expansions  constructed  using  ab-initio  calculations  offer  a 
generalized  potential  that  can  be  used  to  describe  energetics,  since  energies  can  be 
conceived  as  summations  of  the  small  cluster  contributions.  Furthermore,  MBE  technique 
focus  on  positional  degrees  of  freedom,  thus,  it  would  eliminate  a  significant  amount  of 
expensive  and  time  consuming  energy  minimization  required  to  search  for  stable  phase 
structures.  However,  in  practice,  obtaining  the  iV-body  (N> 2)  potentials  is  quite  a 
challenging  problem  and  this  has  been  the  focus  of  our  work. 


Introduction  and  Background 


Advances  in  Molecular  Dynamics  (MD)  and  Monte  Carlo  (MC)  techniques  have  made 
possible  in  the  recent  years  the  systematic  probing  of  material  properties  (phase 
transitions,  thermophysical  properties)  using  computer  experiments.  The  coupling  of 
these  techniques  with  advanced  statistical  methods  would  enable  us  to  systematically 
scan  for  materials  with  extreme  properties.  One  can  imagine  a  scenario  in  which  the 
desired  properties  are  first  specified  and  then  an  extensive  computational  search  is 
performed  to  discover  a  particular  material  that  realizes  those  (Materials  by  Design). 
Subsequently,  targeted  experiments  are  performed  to  actually  create  this  material  in  the 
lab.  Such  a  procedure  would  increase  the  rate  of  new  discoveries  having  a  profound 
effect  on  the  development  of  new  technologies,  saving  at  the  same  time  billions  of  dollars. 
However,  the  magic  ingredient  that  connects  the  MD  and  MC  methods  with  the  real 
world  is  a  physically  accurate  description  of  the  energetics  of  the  materials.  Such  a 
description  is  provided  through  quantum  mechanical  calculations,  albeit  at  a  tremendous 
computational  cost.  Thus,  the  development  of  efficient  ways  to  tabulate  the  results  of 
quantum  mechanical  calculations  is  a  necessary  requirement  towards  the  realization  of 
Materials  by  Design. 

We  first  review  the  atomistic  aspects  of  this  work  using  the  MBE  (multibody  energy 
expansion).  We  want  to  predict  extremal  properties  from  first  principles.  To  do  this 
intelligent  alloying,  we  need  a  method  that  allows  structure  and  property  prediction  of 
multi-atom  systems  that  are  not  necessarily  on  a  Bravais  lattice  (like  FCC,  BCC,  etc.). 
We  need  to  be  able  to  examine  configurations  of  atoms  that  are  placed  anywhere  in  space 
(this  makes  this  method  very  different  from  the  cluster  expansion  method  where  the 
atoms  are  in  a  given  fixed  lattice  and  the  only  thing  you  change  is  what  atom  let  us  say  A 
or  B  you  place  in  each  lattice  location).  We  expand  the  energy  in  two  body,  three -body, 
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etc.  interactions  (Figure  1).  This  is  not  an  obvious  trivial  expansion  as  the  curse  of 
dimensionality  hits  you  fast  and  the  recursive  calculation  of  these  many  body  potentials  is 
very  difficult. 


Approacn 

Require  first-principles  based  calculations.  However,  currently 
infeasible  to  analyze  large  (~104)  systems 

Multi-body  expansions:  Energy'  is  represented  using  a  hierarchy  of 
structure-  independent,  transferable  manv-bodv  potentials 


Construct  these  potentials  from  first-pnnciple  calculations  (VASP, 
Quantum  Espresso  and  DFT++). 


These  potentials  describe  energies  of  arbitrary'  atomic  structures  (~10*) 
as  a  function  of  atomic  positions  very  efficiently 

Figure  1:  The  multibody  energy  expansion  (MBE). 


Methods  for  Constructing  the  Potentials 

In  the  initial  execution  of  this  project,  iV-body  potentials  were  generated  by  tessellating 
the  hyper-surface  and  approximating  the  energy  using  the  finite  element  method. 
Convergence  characteristics  were  significantly  improved  by  weighting  the  energies 
obtained  from  various  truncation  of  the  many  body  expansion.  However,  the  finite 
element  based  tessellation  of  the  hyper-space  places  extensive  restrictions  on  the 
accuracy  of  the  approximation.  Moreover,  the  /V-th  order  potential  lies  in  3N-6 
dimensional  space  and  finite  element  tessellation  (and  subsequent  searching  and 
interpolation)  of  spaces  beyond  6  dimensions  becomes  computationally  ineffective.  In  the 
second  year,  we  incorporated  the  newly  developed  adaptive  sparse  grid  collocation 
(ASGC)  method  based  on  Smolyak  algorithm  into  sampling  the  topology  of  the  clusters 
to  construct  these  iV-body  potentials.  Unfortunately,  we  were  unable  to  interpolate 
energies  of  larger  that  N=5  clusters  with  this  method  because  the  increasing 
dimensionality  of  the  configuration  space  required  a  computationally  forbidding  number 
of  electronic  energy  calculations. 

During  the  third  year,  it  became  apparent  we  had  to  move  to  a  grid-less  interpolation 
scheme.  We  studied  Distance  Geometry  techniques  used  in  the  Protein  Folding  literature 
to  mathematically  describe  the  configuration  space  and  developed  a  new  that  enabled  us 
to  sample  it  efficiently.  The  recently  developed  Multinomial  Expansion  Method  (MEM)  - 
used  in  the  computational  Chemistry  literature  for  the  construction  of  Potential  Energy 
Surface  (PES)  to  study  chemical  reactions  -  was  chosen  as  the  most  promising  candidate 
interpolation  scheme.  It  is  the  first  interpolation  scheme  that  effectively  incorporates  all 
invariance  principles  of  a  potential  energy  surface  in  a  single  functional  form.  The  most 
important  such  invariance  principle  is  the  permutation  invariance  with  respect  to  atoms  of 
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the  same  type.  The  effect  of  this  is  a  drastic  reduction  on  the  number  of  required  ab  initio 
calculations,  thus  making  possible  the  construction  of  MBE  of  higher  orders.  We 
subsequently  improved  the  fitting  capabilities  (MEM)  by  putting  it  in  a  Bayesian 
framework.  The  most  important  new  contribution  is  the  use  of  the  Bayesian  variance  to 
quantify  the  informational  content  of  each  point  in  the  configuration  space  that  lead  us  to 
an  efficient  adaptive  scheme  that  minimizes  the  required  number  of  electronic 
calculations  even  further.  We  demonstrated  that  this  new  technique  considerably 
improves  the  quality  of  the  samples  and  outperforms  the  random  selection  of  data  points. 
We  were  able  to  construct  the  ab  initio  PES  of  Platinum  clusters  of  up  to  6  atoms  with 
only  a  few  thousand  electronic  calculations.  The  ab  initio  PES  were  used  to  find  the 
stable  structures  of  small  Pt  clusters  using  Simulated  Annealing.  The  results  were  found 
to  be  in  very  good  agreement  with  those  found  in  the  literature.  The  constructed  Pt  PES’s 
were  also  used  to  fit  the  interatomic  potentials  up  to  order  6.  It  was  shown  that  those 
become  less  and  less  important  as  their  order  increases,  albeit  slowly  in  low  energy 
regions.  We  used  the  potentials  to  investigate  the  performance  Multi-Body  Expansions  of 
various  orders  for  Pt  clusters  of  up  to  10  atoms.  It  was  demonstrated  that  interactions  of 
at  least  5  atoms  are  required  to  qualitatively  describe  Pt  clusters.  Finally,  we  observed 
that  the  error  introduced  during  the  fitting  procedure  of  the  interatomic  potentials 
propagates  in  a  complicated  manner  through  the  Multi-Body  Expansion  formula  making 
its  naive  application  to  big  clusters  questionable.  It  is  the  object  of  our  current  research  to 
investigate  the  propagation  of  this  error  through  the  MBE  formula  and  design  effective 
techniques  to  filter  it  out.  We  believe  that  such  filtering  schemes  have  to  be  case  specific 
(different  for  each  material)  and  should  utilize  further  physical  information.  This  problem 
constitutes  the  final  obstacle  towards  the  construction  of  fully  transferable  potential 
energy  surfaces  using  the  MBE  framework. 

With  the  PES  surfaces  in  place,  exploration  of  the  energy  landscape  in  the  high 
dimensional  configuration  space  becomes  an  easier  task  that  could  potentially 
revolutionize  the  search  for  materials  for  extremal  properties.  In  our  immediate  plans  we 
are  working  towards  integration  of  PES  meta  models  with  MD  and  MC  techniques  to 
allow  us  computing  materials  with  desired  mechanical  and  thermophysical  properties, 
phase  transition  characteristics,  etc.  Many  applications  to  surface  design  (e.g.  of  Pt 
clusters  to  maximize  H-adsorption)  are  also  anticipated. 

In  the  remaining  of  this  final  report,  we  first  briefly  discuss  the  developments  of  FEM 
like  tessellation  techniques  of  the  configuration  space  for  PES  construction.  For  brevity 
of  the  report  we  do  not  discuss  the  activities  on  the  sparse  grid  interpolation  approach 
since  the  number  of  ab  initio  simulations  needed  with  this  method  was  prohibitory  high. 
We  finally  conclude  with  the  Bayesian  framework  for  interpolating  potentials  using 
invariant  polynomials  in  the  high-dimensional  configuration  space.  It  is  this  framework 
that  we  believe  provides  the  best  available  option  for  PES  surrogate  construction  using 
the  minimum  number  of  ab  initio  runs  for  properly  chosen  realizations  in  the 
configuration  space. 
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The  effect  of  structural  relaxations  in  alloys  is  described  using  a  multi-body  energy  expansion 
formalism.  IV-body  potentials  in  the  multi-body  expansion  are  computed  from  energies  of  isolated 
clusters,  which  in  turn,  are  calculated  from  empirical  potentials  or  self-consistent  quantum  mechan¬ 
ical  calculations.  Convergence  characteristics  of  multi-body  expansions  (MBE)  are  improved  by 
weighting  energies  obtained  from  various  truncations  of  many-body  expansion  in  a  new  method 
called  weighted  MBE  (wMBE).  It  is  shown  that  multi-body  expansions  of  many-atom  systems  can 
be  efficiently  constructed  using  interpolation  of  isolated  cluster  energies  from  databases.  In  contrast 
to  the  method  of  cluster  expansion,  wMBE  focuses  on  positional  degrees  of  freedom  and  hence, 
explicitly  handles  structural  relaxations  during  computations  of  stable  atom  clusters  and  periodic 
or  amorphous  phase  structures. 
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I.  INTRODUCTION 

Calculation  of  stable  structures  of  alloys,  clusters,  sur¬ 
faces  and  molecules  from  first-principles  is  an  important 
step  towards  design  of  materials  with  exceptional  proper¬ 
ties.  Identification  of  stable  alloy  phases  aid  in  construc¬ 
tion  of  phase  diagrams  from  first-principles.  Because  of 
the  immense  variety  of  phase  structures,  identification  of 
stable  structures  at  different  combinations  of  the  alloying 
elements  is  a  non-trivial  problem.  While  a  first-principles 
approach  based  on  density  functional  theory  (DFT)  pro¬ 
vide  a  rigorous  way  for  calculating  formation  energies  of 
phase  structures,  the  computational  complexity  of  per¬ 
forming  fully-relaxed  calculations  over  the  entire  set  of 
possible  phase  structures  makes  this  method  prohibitive. 
Techniques  such  as  cluster  expansion  1-7  and  more  re¬ 
cently,  data  mining  techniques8,9  allow  one  to  accelerate 
the  search  for  stable  phase  structures. 

In  cluster  expansion  methods  (CEM)1-4,  the  relaxed 
energy  of  an  atomic  structure  is  represented  as  a  lin¬ 
ear  combination  of  characteristic  energies  of  clusters  of 
atoms  over  a  fixed  lattice.  The  coefficients  in  the  clus¬ 
ter  expansion  are  computed  using  relaxed  DFT  energy 
calculations  of  few  prototype  structures1.  This  method 
includes  only  ordering  degrees  of  freedom  as  provided  by 
different  possible  arrangements  of  atom  types  on  a  fixed 
parent  lattice.  Consequently,  CEM  fails  in  cases  where 
the  alloy  phases  have  complex  structures  that  are  dif¬ 
ferent  from  the  superstructures  of  the  underlying  parent 
lattice  (for  example,  FCC  or  BCC  lattices)  and  exhibits 
convergence  issues  in  cases  where  structural  relaxation 
effects  are  dominant5,6  (for  example,  in  alloys  involving 
constituents  with  large  size  differences). 

In  another  technique  called  multi-body  expansion 


(MBE),  TV-body  potentials  (or  otherwise,  cluster  poten¬ 
tials10)  constructed  from  ab-initio  calculations  are  used 
to  describe  energies  of  arbitrary  atomic  structures  as  a 
function  of  atom  positions.  The  total  energy  is  repre¬ 
sented  as  a  summation  over  potentials  of  underlying  iso¬ 
lated  atom  clusters  in  the  structure,  with  series  terms  in¬ 
volving  pair,  three-body,  four-body,..,  IV-body  potentials. 
Up  to  third-order  truncations  of  multi-body  expansions 
have  been  previously  used  in  related  empirically  derived 
potentials,  namely  the  Gupta11  and  Murrell-Mottram 
(MM) 12-15  potentials.  Multi-body  potentials  focus  on 
positional  degrees  of  freedom  and  hence,  explicitly  han¬ 
dles  structural  relaxations  during  computations  of  sta¬ 
ble  phase  structures.  Structural  relaxation  effects  can 
be  treated  in  a  cluster  expansion  approach  by  combin¬ 
ing  it  with  position-dependent  potentials  in  the  form  of 
a  hybrid  cluster  expansion6.  Another  method  combining 
CEM  and  multi-body  potentials  was  proposed  recently 
for  introducing  positional  degrees  of  freedom  in  a  more 
generalized  cluster  expansion16,17.  However,  building  the 
IV-body  potentials  from  atomistic  calculations  is  quite 
a  challenging  problem.  Firstly,  the  number  of  clusters 
(and  thus,  the  number  of  cluster  energies  that  need  to  be 
calculated)  in  an  IV-atom  system  increases  geometrically 
with  the  order  of  expansion  (Fig.  1).  Secondly,  although 
the  approach  provides  good  convergence  for  rare-gas  crys¬ 
tals,  convergence  of  MBE  is  not  smooth  for  metallic  crys¬ 
tals18,21.  Absence  of  smooth  convergence  does  not  allow 
establishment  of  a  hard  cut-off  for  the  series  terms.  Con¬ 
sequently,  there  has  been  no  published  reports  of  a  multi¬ 
body  expansion  constructed  directly  from  first-principle 
calculations. 

This  paper  addresses  these  drawbacks  by  proposing  a 
multi-body  expansion  with  weighted  terms.  Using  this 
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FIG.  1:  Increase  in  the  total  number  of  clusters  involved  in  a 
multi-body  expansion  for  a  32-atom  system  with  the  order  of 
expansion. 


technique,  it  becomes  possible  to  accurately  compute  en¬ 
ergies  of  iV-atom  systems  from  knowledge  of  small  cluster 
energies  computed  from  first  principles.  The  efficiency 
of  this  new  method  called  weighted  multi-body  expan¬ 
sion  (wMBE)  is  emphasized  through  examples  in  this 
work.  Another  contribution  of  this  paper  is  a  formal 
technique  to  rapidly  calculate  multi-body  expansions  us¬ 
ing  linear  interpolation  over  tessellated  cluster  configura¬ 
tional  spaces.  MBE  using  interpolated  energies  is  several 
orders  of  magnitude  faster  than  using  DFT  calculations, 
since  cluster  energies  are  computed  beforehand  and  are 
directly  sampled  from  the  database  when  computing  the 
multi-body  expansion. 


A.  Multi-body  expansion  methodology 


Consider,  for  instance,  a  configuration  of  M- atoms 
(possibly  all  different),  whose  energy  we  intend  to  com¬ 
pute.  We  denote  the  total  energy  of  this  M-particle 
system  using  Ep  =  Ep{X\,  X2,  ■ .  . , !«),  where  P  is 
the  order  of  the  expansion  used  and  the  position  Rn 
of  atom  n  is  grouped  together  with  the  species  of  atom 
n  denoted  by  an  integer  an ,  Xn  =  (R„,  <rn).  As  the 
order  of  labelling  the  M  atoms  is  arbitrary,  the  form 
of  Ep(X i, . . . ,  Xi, . . . ,  Xj, . . . ,  Xm )  must  be  symmetric 
with  respect  to  interchange  of  Xi  and  Xj . 

From  here  on,  we  denote  M  as  the  total  number  of 
atoms  in  the  system,  N  =  1,2,  ...,P  denotes  a  N- 
atom  cluster  within  the  M— atom  system.  Further,  L  = 
1,  2, . . . ,  N  denotes  an  arbitrary  L-atom  cluster  within  an 
A— atom  cluster.  The  energy  Ep  of  an  M-particle  sys¬ 
tem  is  represented  as  summation  over  a  series  of  Al-body 


interaction  potentials  V^N'>  via 

p 

Ep(X1,X2,...,Xm)  =  ^2eW(X1,X2,...,Xm), 

N=  1 

MM  M 

^W=E  E 

mi  =  lm2=mi  +  l  mjv=mjv- i  +  l 

v^N\xmi,xm2,...,xmN).  (1) 

The  potentials  can  be  inverted  via  the  Mobius  inver¬ 
sion  approach  from  number  theory.  Mobius  inversion  has 
been  used  previously  for  extraction  of  potentials  from  en¬ 
ergy  data  by  Chen19,20  although  in  a  different  context. 
In  the  case  of  multi-body  potentials  V^N\  the  Mobius 
inversion  is  given  as16: 


N 


N 


yW(x1,x2,...,x]V)  =  ^(-i)JV-i  Y, 

L— 1  mi— 1 

N  N 

E  •••  E  E*(Xmi,Xm2,...,XmL)(  2) 


m2=mi  +  l  mL—rriL- i  +  l 


Here,  we  denote  the  energies  of  L— atom  clusters  within 
the  A— atom  cluster  as  E* .  The  above  equation  con¬ 
stitutes  a  unique  definition  of  iV-body  potentials 
which  are  structure-independent  because  this  equation 
does  not  carry  any  information  about  the  environment 
of  the  atom  clusters16.  V^2\Xi,Xj)  can  be  understood 
as  the  excess  energy  attributed  to  pair  interactions  in  an 
isolated  atom  pair  i,j,  i.e.,  V^2\Xi,Xj)  =  E*(Xi,Xj)  — 
E*  ( Xj )  —  E*  (Xj ) .  Similarly,  (Xi ,  Xj ,  Xu )  can  be  un¬ 
derstood  as  the  excess  energy  attributed  to  three-body 
interactions  in  a  isolated  trimer  (i,j,  k ): 

V^(Xi,  Xj,  Xk)  =  E*(Xi,Xj,Xk) 

-  (VM(Xi,Xj)  +  vW(Xj,Xk)  +  V^(Xi,Xk)) 

-  (■ E*(Xi )  +  E*(Xj)  +  E*(Xk)).  (3) 

Once  the  potentials  have  been  constructed,  they 
can  be  used  to  calculate  the  energy  Ep(X \,X2,  ■  ■  ■ ,  Xm) 
for  a  M-atom  system  using  Eq.  (1).  The  first  critical 
requirement  of  the  technique  is  the  knowledge  of  com¬ 
plete  energy  surface  (cluster  energies  versus  atom  po¬ 
sitions  and  types)  of  small  isolated  clusters  of  atoms 
(E* ,  L  =  1, . . . ,  5)  for  building  the  potentials  in  Eq.  (2). 
Secondly,  it  is  essential  that  the  expansion  converges 
within  a  small-order  of  expansion  (i.e.  P  <  5)  for  com¬ 
putational  efficiency.  These  two  aspects  are  addressed  in 
the  next  two  sections.  Complete  energy  surface  for  small 
isolated  clusters  is  created  by  mathematically  defining 
the  configurational  space  of  clusters,  tessellation  of  the 
space,  computation  of  cluster  energies  on  nodal  points, 
followed  by  interpolation  of  cluster  energies  as  described 
in  the  next  section.  Computational  efficiency  is  improved 
through  the  use  of  weighted  multi-body  expansions  as  ex¬ 
plained  in  section  C.  Efficiency  can  be  further  improved 
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by  performing  computations  in  parallel  by  distributing 
the  M  atoms  involved  in  the  loop  index  mi  in  Eq.  (1)  to 
different  processors. 


B.  Construction  of  cluster  energy  surfaces 

The  basic  idea  of  the  approach  to  rapidly  compute 
multi-body  expansions  of  arbitrary  systems  is  to  build 
an  interpolation  function  for  the  isolated  cluster  ener¬ 
gies  E*  from  the  pre-computed  database.  Given  a  set 
of  n  ?n-atom  clusters  represented  as  0  =  {Cd}i=i  in  the 
d-dimensional  configurational  space,  we  try  to  build  a 
smooth  function  that  maps  clusters  to  ab-initio  energies, 
/  :  Rd  — >  K.  In  particular,  we  use  an  interpolant  If  such 
that  !/(&)=/(£,),  V*=l,...,n. 

The  first  step  in  this  procedure  will  be  to  define  the 
d-dimensional  configurational  space  of  an  ?n-atom  clus¬ 
ter.  The  positions  of  the  atoms  in  the  cluster  are  rep¬ 
resented  by  the  distance  between  atoms,  Rij  >  0.  For 
two-atom  clusters  (m  =  2),  the  configurational  space 
is  one-dimensional,  with  each  point  x  in  the  space  rep¬ 
resenting  a  two-atom  cluster  with  inter-atomic  distance 
of  I?i2  =  x.  As  the  number  of  atoms,  ?n,  in  the  clus¬ 
ter  increases,  the  number  of  distances,  Rij  necessary  to 
completely  and  uniquely  describe  the  cluster  increases 
rapidly.  Up  to  m  <  4,  clusters  are  uniquely  represented 
by  ^m(m  —  1)  independent  variables. 

For  example,  the  space  of  all  possible  three-atom  clus¬ 
ters  is  three-dimensional  as  shown  in  Fig.  2(a).  This 
space  is  a  convex  hull  with  9  planes  (symmetries  not  in¬ 
cluded)  due  to  a  linear  set  of  constraints  arising  from 
three  triangle  inequalities  of  the  form  R,tj  +  Rjk  >  Rik 
that  constrain  the  location  of  atoms  in  the  three-atom 
cluster  and  the  upper  and  lower  cutoff  used  for  possible 
cluster  sizes  in  the  database:  Rij  >  l  and  Rij  <  u  with 
i,  j,  k  =  1 ...  3.  Cluster  symmetries  can  be  used  to  further 
reduce  the  space  and  consequently,  reduce  the  number  of 
energy  calculations  required.  Figure  2(b)  shows  the  re¬ 
duced  space  accounting  for  symmetries  (R12  <  R23  < 
R13)  in  the  case  where  all  3  atoms  are  of  the  same  type. 
Also  shown  in  Fig.  2(a)  is  the  tessellation  of  the  config¬ 
urational  space  of  clusters.  The  energy  of  a  cluster  cor¬ 
responding  to  each  nodal  point  in  the  space  is  calculated 
and  stored  in  the  database.  The  plot  of  energy  versus  in¬ 
teratomic  distance  for  a  two-atom  Pt  cluster  is  shown  in 
Fig.  3  with  location  of  nodal  points  for  two-atom  clusters. 
Higher-dimensional  spaces  are  adaptively  tessellated  as 
shown  in  Fig.  2(a)  with  a  finer  discretization  of  regions  in¬ 
volving  small  clusters.  The  tessellation  of  the  configura¬ 
tional  space  is  carried  out  using  n-dimensional  Delaunay 
triangulation,  as  implemented  in  the  qhull22  program. 
Tessellation  generates  elements  (known  as  a  simplex)  over 
which  local  linear  interpolation  is  carried  out  to  find  the 
energy  of  any  other  three-atom  cluster  within  the  space. 
The  discretization  and  interpolation  techniques  are,  in 
essence,  same  as  those  used  in  the  popular  finite  element 
techniques  for  PDEs.  Energy  (E)  of  an  arbitrary  cluster 


FIG.  2:  (left)  shows  the  space  of  all  possible  three-atom  clus¬ 
ters  within  an  upper  and  lower  cutoff  cluster  size.  This  space 
represents  a  convex  hull  in  3D.  (right)  Use  of  symmetries  (in 
the  case  where  all  three  atoms  are  of  one  type,  e.g.  Pt-Pt-Pt 
clusters)  can  further  reduce  the  space.  The  simplices  used  to 
perform  local  linear  interpolation  of  energies  are  also  shown. 
In  3D,  the  simplex  is  a  tetrahedron. 
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FIG.  3:  The  plot  of  energy  versus  interatomic  distance  for  a 
two-atom  Pt  cluster.  The  location  of  nodal  points  in  this  one¬ 
dimensional  configurational  space  and  the  upper  cutoff  used 
for  calculations  are  indicated. 


with  cluster  specifier  =  [£*,£2,  ••  • , Q]  in  a  tessellated 
d-dimensional  configurational  space  is  given  as: 

E  =  aT  Ee,  (4) 

where  Ee  is  the  vector  containing  energies  at  the  nodes 
of  the  simplex  within  which  £2  is  located,  a  is  ob¬ 
tained  as  a  =  A  xb,  where  A  =  [1,  . . . ,  ]T  and 

b  =  [1,  £1 ,  £2 !  ■  •  •  i  Q]T ■  Here,  ^  denotes  a  vector  contain¬ 
ing  the  ith  coordinate  value  of  all  nodes  in  an  element 
e.  The  element  e  is  located  by  calculating  a  for  every 
element  in  sequence  and  selecting  the  element  e  where 
all  elements  of  a  >  0.  This  step  becomes  more  time- 
consuming  as  the  dimensionality  of  configurational  space 
(hence,  the  number  of  elements)  increases.  Further,  the 
geometry  of  the  configurational  space  becomes  more  com¬ 
plex  as  the  dimensionality  of  the  configurational  space 
increases.  For  example,  the  configurational  space  of  a 
fourth-order  cluster  (excluding  symmetries)  involves  24 
linear  constraints  and  a  quartic  constraint. 

The  number  of  independent  variables  specifying  a  m  > 
4  atom  cluster  is  given  by  d  =  3m  —  6  although  )m(m- 1) 
variables  are  needed  to  uniquely  define  a  cluster23.  An 
example  of  how  cluster  specifiers  are  determined  for  a 
5-atom  cluster  is  illustrated  in  Fig.  4.  In  this  example, 
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FIG.  4:  In  the  case  of  a  five-atom  cluster,  the  locations  of 
the  fourth  and  fifth  atoms  can  be  fixed  with  respect  to  the 
plane  formed  by  atoms  1  —  2  —  3  using  the  cluster  specifiers 
[f?i2,  R23,  R13,  R14,,  R24,  R34,  R15,  R25,  -R35].  However,  these 
specifiers  do  not  completely  represent  the  cluster.  The  de¬ 
pendent  variable  in  this  case  is  R45  which  can  take  one  of 
possible  two  values  based  on  the  location  of  atom  5  either 
above  or  below  the  plane  formed  by  atoms  1  —  2  —  3. 


there  are  9  independent  variables  and  1  dependent  vari¬ 
able  (R45)  that  can  take  one  of  two  values  based  on  the 
location  of  the  fifth  atom.  Thus,  m  >  4  cases  present 
special  difficulties  associated  with  dependent  variables. 
We  address  this  issue  by  creating  different  configuration 
spaces  corresponding  to  the  values  that  each  dependent 
variable  takes.  For  the  case  of  a  5-atom  cluster,  this 
means  that  two  potentials  need  to  be  created,  one  for  the 
case  where  atom  5  is  above  the  plane  formed  by  atoms 
1  —  2  —  3  and  another  when  it  is  below  that  plane. 

For  a  binary  AB  system,  all  possible  cluster  con¬ 
figurational  spaces  are  created  for  a  given  clus¬ 
ter  size,  e.g.  for  L-atom  clusters,  L  +  1  energy 
databases  (e.g.  for  L  =  2,  3  databases  containing 
E*(Xa,  Xa),  E*(Xa,  Xb),  E*(Xb,  Xb))  need  to  be  gen¬ 
erated.  The  upper  and  lower  cutoff  were  selected  by 
carefully  analyzing  the  energies  of  two-atom  clusters  over 
a  large  range  of  B12  to  locate  an  upper  cut-off  beyond 
which  the  interaction  between  atoms  were  not  significant 
and  a  lower  cutoff  where  the  interaction  energy  is  posi¬ 
tive. 

For  Platinum  with  lattice  parameter  of  a  =  7.5  bohr, 
the  lower  cutoff  of  atom  spacing  in  a  cluster  within  the 
database  was  fixed  as  Rij  >  0.3a  and  upper  cutoff  was 
fixed  as  Rij  <  1.5a.  The  plot  of  energy  versus  inter¬ 
atomic  distance  for  a  two-atom  Pt  cluster,  from  which 
the  cut-offs  were  identified,  is  shown  in  Fig.  3.  The  cut¬ 
offs  signify  that  clusters  with  R,j  <  0.3a  and  R,VI  >  1.5a 
are  not  available  in  the  database.  During  MBE  calcu¬ 
lations,  energies  of  clusters  containing  such  interatomic 
distances  are  approximated  using  the  following  means. 
For  R^  <  0.3a,  cluster  energies  were  given  an  artificial 
high  value  to  signify  that  such  configurations  are  not  en¬ 
ergetically  feasible.  For  Af-atom  clusters  with  R^  >  1.5 a, 
the  excess  energy  attributed  to  N-body  interactions  is  as¬ 
sumed  to  be  zero  (i.e.  V ^  ~  0,  for  iV-atom  clusters  with 
an  Rij  >  1.5a).  This  is  mathematically  equivalent  to  ap¬ 
proximating  the  energies  of  large  clusters  using  energies 
of  smaller  sub-clusters.  For  example,  Fig.  5(a)  shows  the 
energy  surface  of  three-atom  Platinum  clusters  up  to  an 


(a)  (b) 

FIG.  5:  Energy  surface  (E*(X  1,  X2,  X3))  for  3-atom  Pt  clus¬ 
ters  whose  atoms  are  positioned  at  the  vertices  of  a  right  an¬ 
gled  triangle  with  line  joining  atoms  A'2  and  A3  forming  the 
hypotenuse.  Figure  (a)  shows  computed  Platinum  three-atom 
cluster  energies,  while  (b)  shows  extension  of  energies  beyond 
the  cutoff  using  energies  of  smaller  clusters  (E*( Xi,Xj)  and 

E*(  A'O). 

upper  cutoff  size  of  1.5a  and  Fig.  5(b)  shows  the  complete 
energy  surface  when  the  energies  beyond  the  upper  cutoff 
are  approximated  using  two-  and  one-atom  energies. 

C.  Weighted  multi-body  expansion 

Multi-body  expansion  has  been  shown  to  work  very 
well  for  rare-gases  where  the  expansion  is  dominated  by 
pair  interactions  making  higher  terms  in  the  expansion 
negligible.  Total  energy  of  metallic  systems,  however, 
has  significant  contributions  from  higher-order  interac¬ 
tions  and  the  expansion  has  non-smooth  convergence  be¬ 
havior  18 .  Figure  6  shows  the  behavior  of  multi-body 
expansion  for  an  eight  atom  (2  unit  cell)  FCC  Platinum 
cluster  that  requires  at  least  a  7th  order  expansion  to 
capture  the  true  energy.  It  is  observed  here  that  ener¬ 
gies  computed  by  including  successively  higher-orders  of 
interaction,  in  fact,  oscillate  around  the  true  energy.  An 
ad-hoc  numerical  approach  for  estimating  the  true  ener¬ 
gies  for  the  case  in  Fig.  6  will  be  to  appropriately  weight 
the  energies  obtained  at  different  orders  of  multi-body 
expansion,  which  is  akin  to  smoothing  (or  filtering)  the 
energy  oscillations  in  Fig.  6.  Numerical  experiments  pre¬ 
sented  in  the  next  section  indicate  that  weighted  MBE 
(wMBE)  calculations  lead  to  dramatic  improvement  to 
the  convergence  behavior  of  the  multi-body  expansion. 
In  the  wMBE  approach,  the  energies  up  to  a  cut-off  or¬ 
der  of  expansion  P  are  weighted  so  that  we  reach  as  close 
to  the  true  energy  (Em)  of  an  M- atom  system  as  follows: 

Em(X  1,  A'2, . .  • ,  Xm)  =  a\Ei(Xi,X2, ,  Xm) 

+  a2E2(Xi , . . . ,  Xm)  +  ■  ■  ■  +  apEp(X  1, . . . ,  A'm)(5) 

The  coefficients  a  =  [aq,  a2,  ■  ■  . ,  ap]T  are  computed 
by  solving  the  equation: 

a  =  C+E ,  (6) 

where,  E  are  the  true  energies  of  q  M- atom  clusters 
(Xl,i  =  1  ,...,q)  computed  with  self-consistent  DFT 
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calculations.  Each  row  of  C  contains  the  energies 
[E±,  E‘2, ,  Ep ]W  obtained  from  multi-body  expansion 
of  each  of  these  clusters  (where  Ep^'  =  Ep(X1)).  C+  is 
the  pseudo-inverse  of  matrix  C.  The  technique  to  obtain 
coefficients  a  is  thus,  similar  to  the  method  of  Connolly 
and  Williams1  used  for  cluster  expansions.  In  their  tech¬ 
nique,  truncation  of  the  expansion  is  based  on  which  clus¬ 
ters  are  important,  for  example,  in  FCC  crystals  where 
only  clusters  containing  nearest  neighbors  are  important, 
the  series  is  truncated  at  fourth-order.  In  the  case  of  the 
wMBE,  however,  the  cutoff  of  the  order  of  interactions 
needs  to  be  identified  through  numerical  experiments  as 
will  be  demonstrated  in  the  next  section. 
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FIG.  6:  Convergence  of  the  many-body  energy  expansion  of 
an  eight-atom  FCC  Platinum  cluster  requires  at  least  a  7th 
order  expansion  to  reasonably  capture  the  true  energy. 


Before  proceeding  to  examples,  we  summarize  the  steps 
involved  in  the  overall  algorithm  as  follows: 

1.  Offline  calculations  (steps  1-4):  Constructing  a 
database. 

Generate  coordinates  ({£d)T=i)  for  sampling  the 
configurational  space  of  all  cluster  sizes  involved. 
For  example,  each  three-atom  cluster  (1-2-3)  cor¬ 
responds  to  the  coordinate  Q  =  (R12,  R23,  R13)  in 
the  configurational  space  where  e.g.  R12  is  the  in¬ 
teratomic  distance  between  atoms  1  and  2.  During 
this  step,  various  constraints  based  on  geometry  or 
symmetry  are  used  to  reduce  the  number  of  nodes 
in  the  configurational  space. 

2.  For  an  L— atom  cluster  of  a  binary  system,  coordi¬ 
nates  for  L  +  1  configurational  spaces  need  to  be 
created  during  step  1.  Each  configurational  space 
corresponds  to  a  different  atom  type  list,  e.g.  for 
L  =  2  atom  clusters  of  a  binary  alloy  A  —  B,  config¬ 
urational  spaces  for  clusters  of  types  A  —  A,  A  —  B 
and  B  —  B  need  to  be  generated. 

3.  Perform  tessellation  of  coordinates  in  all  config¬ 
urational  spaces  and  store  nodal  coordinates  and 
element-node  fists  in  the  database. 

4.  Generate  input  files  and  perform  self-consistent 
DFT  calculations  to  compute  energies  {E*(£ld))  at 


nodal  locations  of  all  configurational  spaces.  Ener¬ 
gies  from  the  DFT  calculation  are  read  and  stored 
in  databases,  one  corresponding  to  each  configura¬ 
tional  space. 

5.  Calculation  of  MBE  coefficients. 

Compute  self-consistent  ab-initio  calculations  to 
compute  energies  (to  obtain  E  in  Eq.  (6))  of  a  few 
(three  or  four)  different  IV-atom  configurations. 

6.  Compute  energies  using  MBE  with  increasing  or¬ 
ders  of  expansion  and  obtain  [Ei,  E2, ... ,  Ep ]W  for 
each  IV-atom  configuration  used  in  step  5.  During 
multi-body  expansion,  potentials  (T/W)  are  cre¬ 
ated  using  Eq.  (2)  on  the  fly,  using  cluster  energies 
E*  obtained  by  interpolating  from  the  database 
constructed  in  steps  (1-4).  The  steps  involved  to 
compute  cluster  energy,  E* ,  of  an  arbitrary  cluster 
are: 

(a)  Locate  the  cluster  in  the  corresponding  config¬ 
urational  space.  For  example,  a  three-atom  cluster 
of  type  A  —  A  —  B,  is  located  at  the  coordinate 
(f  =  {R\2,  R23,  R13)  in  the  three  dimensional  con¬ 
figurational  space  of  A  —  A  —  B  type. 

(b)  Identify  the  element  in  which  the  cluster  is  lo¬ 
cated  and  perform  linear  interpolation  using  known 
energies  at  nodal  values  in  that  element  using 
Eq.  (4). 

(c)  Energies  of  clusters  that  are  not  available  in 
the  database  are  approximated  using  the  methods 
detailed  at  the  end  of  Section  B. 

7.  Compute  the  coefficients  a  of  weighted  MBE  using 
Eq.  (6).  Perform  tests  for  convergence  by  compar¬ 
ing  energies  predicted  by  wMBE  with  ab-initio  cal¬ 
culation  for  few  other  configurations  of  iV-atoms. 

8.  MBE  calculations  for  arbitrary  N-atom  configura¬ 
tions. 

The  converged  weighted  expansion  can  be  now  be 
employed  for  computing  energies  of  other  iV-atom 
systems  using  Eq.  (5)  and  Eq.  (1).  During  calcu¬ 
lations,  cluster  energies  E*  are  again  interpolated 
from  the  database  as  in  step  6. 


D.  Results  for  metallic  systems 

1.  Extrapolatory  performance  of  wMBE  approach :  In 
the  first  test  case,  energies  predicted  by  multi-body  ex¬ 
pansion  are  compared  with  true  energies  obtained  using 
the  embedded  atom  potential  of  Sutton  and  Chen24  for 
Platinum  atom  clusters.  Convergence  of  the  expansion  is 
tested  using  exact  cluster  energies  (without  performing 
interpolation).  Atom  configurations  used  in  these  cases 
correspond  to  nx  x  ny  x  nz  clusters  with  rq  unit  cells 
located  in  the  ith  direction. 

Figure  7  shows  the  energies  obtained  for  an  isolated 
4x1x1  (16-atom)  cluster  of  Platinum  computed  using 
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FIG.  7:  Comparison  of  true  energies  of  a  4  x  1  x  1  (16-atom) 
FCC  Platinum  cluster  with  that  predicted  by  multi-body  ex¬ 
pansion.  Weights  in  the  multi-body  expansion  were  computed 
using  3  energies  at  lattice  parameters  of  7.2,  7.4  and  7.6  bohr. 


FIG.  8:  Comparison  of  true  energies  of  a  4  x  1  x  1  (16-atom) 
FCC  Platinum  cluster  with  MBE  expansion  results  for  an 
extrapolatory  case  of  a  FCC  lattice  where  the  face  centered 
atom  in  the  x  —  y  plane  and  y  —  z  plane  in  the  FCC  basis  are 
translated  by  (—0.1, —0.1,  0)  and  (0,0.1,  0.1),  respectively  in 
crystal  coordinates. 


2nd,  3rd  and  4th  order  wMBE  and  the  true  energies.  In 
all  cases,  the  parameters  a.  were  computed  using  4x1x1 
Pt  clusters  using  just  3  energies  at  lattice  parameters  of 
7.2,  7.4  and  7.6  bohr  as  indicated  in  Fig.  7.  Energies  in¬ 
crease  linearly  within  this  range  of  lattice  parameters.  In 
spite  of  this,  predicted  energies  from  the  3rd  and  4th  or¬ 
der  multi-body  expansions  exactly  capture  the  parabolic 
nature  of  the  true  energy  profile.  As  a  test  of  the  ex¬ 
trapolatory  performance  of  wMBE,  we  perturb  the  face 
centered  atoms  in  the  x  —  y  plane  and  y  —  z  plane  of 
the  FCC  basis  by  (—0.1,  —0.1, 0)  and  (0, 0.1, 0.1),  respec¬ 
tively  in  crystal  coordinates.  In  spite  of  the  large  changes 
in  energy  resulting  from  this  perturbation,  the  expansion 
built  previously  for  a  FCC  cluster  is  able  to  reproduce  the 
energy  profile  of  this  distorted  cluster  accurately  (Fig.  8). 

2.  Convergence  of  wMBE  in  extrapolatory  cases :  Fig¬ 
ure  9  shows  the  energies  predicted  at  various  lattice  pa¬ 
rameters  using  2nd  and  3rd  order  wMBE  for  an  isolated 
2x2x1  FCC  Platinum  cluster.  In  this  case,  the  param¬ 
eters  a  were  originally  computed  using  2x2x1  FCC 
clusters  of  Pt  using  11  lattice  parameters  between  6  bohr 
to  8  bohr  in  the  increments  of  0.2  bohr.  Although  Fig.  9 


FIG.  9:  Comparison  of  true  energies  obtained  for  a  2  x  2  x 
1  (16-atom)  FCC  Platinum  cluster  with  energies  computed 
using  2nd  and  3rd  order  wMBE. 


shows  that  the  third-order  expansion  is  adequate  to  cap¬ 
ture  the  true  energy  profile,  use  of  higher-orders  of  ex¬ 
pansion  improve  the  performance  in  extrapolatory  cases. 
Figure  10  depicts  the  performance  of  3rd,  4th  and  5th  or¬ 
der  MBE  in  an  extrapolatory  case  where  the  face  centered 
atom  in  the  x  —  y  plane  of  the  FCC  basis  is  translated  by 
(—0.1, —0.1, 0)  in  crystal  coordinates.  Figure  11  shows 
the  decrease  in  the  l 2  norm  error  in  energies  predicted 
with  increasing  order  of  multi-body  expansion.  Several 
other  numerical  experiments  of  this  kind  indicate  that  the 
wMBE  approach  captures  the  energy  profile  for  any  ran¬ 
dom  configurations  of  iV-atom  Pt  clusters,  and  thus,  has 
potential  applications  in  NVE  or  NVT  atomistic  simula¬ 
tions.  The  weighting  procedure  aims  to  average  out  the 
extraneous  energy  contributions  (eg.  surface  energies) 
arising  due  to  lack  of  environment  in  isolated  clusters. 
The  limitation  in  the  procedure  is  that  a  change  in  num¬ 
ber  of  atoms  (TV)  simulated  necessitates  re-calibration  of 
MBE  coefficients.  Figure  9  shows  the  energy  variation 
with  lattice  parameter  obtained  from  an  MBE  expan¬ 
sion  calculated  using  cluster  energies  interpolated  from 
a  database.  For  interpolation,  the  second-order  configu¬ 
rational  space  is  discretized  into  20  linear  elements  (21 
nodes)  and  the  third-order  configurational  space  (includ¬ 
ing  symmetries)  was  approximated  using  16374  tetrahe¬ 
dral  elements  on  which  energies  were  calculated  at  3191 
nodal  locations.  Although  discretization  and  linear  in¬ 
terpolation  introduce  errors  in  calculation  of  energies,  it 
is  seen  that  the  technique  still  reasonably  captures  the 
energy  profile  and  the  energy  minima.  The  advantage 
of  interpolation  approach  is  that  it  is  order  of  magni¬ 
tude  faster  since  cluster  energies  are  computed  before¬ 
hand  and  are  directly  sampled  from  the  database  during 
simulations. 

3.  wMBE  using  interpolated  energies  from  ab-initio 
calculations :  Figure  12  shows  structure  optimization 
to  find  the  lattice  constants  for  FCC  Platinum  system 
using  interpolated  energies  of  clusters  computed  from 
first  principles  DFT  calculations.  Since  MBE  inher¬ 
ently  uses  non-periodic  configurations,  energies  of  peri¬ 
odic  structures  are  computed  by  considering  supercells 
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FIG.  10:  Comparison  of  true  energies  obtained  for  a  2  x  2  x 
1  (16-atom)  FCC  Platinum  cluster  with  energies  computed 
using  3rd,  4th  and  5th  order  wMBE  for  an  extrapolatory  case 
of  a  FCC  lattice  with  the  face  centered  atom  in  the  x—y  plane 
of  the  FCC  basis  is  translated  by  (—0.1,  —0.1,  0)  in  crystal 
coordinates. 


FIG.  11:  Decrease  in  the  l2  norm  error  in  energies  with  in¬ 
creasing  order  of  multi-body  expansion  for  the  extrapolatory 
case  in  Fig.  10. 


with  true  energies  E  used  for  fitting  Eq.  (5)  obtained 
from  self-consistent  DFT  calculations  of  a  periodic  unit 
cell.  In  this  example,  a  5  x  5  X  5  (500-atom)  FCC 
cluster  is  considered.  The  variation  of  cohesive  energy 
(EC(X =  E*(Xl7..,Xm) -J2ZiE*(Xi))  of  3- 
atom  clusters  with  interatomic  distances  (R12, -R23, -R13) 
is  shown  on  the  configurational  space  (accounting  for 
symmetry)  in  Fig.  12.  The  configurational  space  is  dis¬ 
cretized  into  4609  tetrahedral  elements  on  which  linear 
interpolation  is  carried  out.  Ab-initio  energy  data  were 
computed  on  1027  nodal  locations.  Figure  13  shows  com¬ 
parison  of  the  energies  computed  using  3rd  and  4th  order 
wMBE  with  the  true  energies.  Coefficients  in  the  multi¬ 
body  expansion  were  generated  using  three  ab-initio  en¬ 
ergy  calculations  of  a  periodic  FCC  Platinum  lattice  with 
lattice  parameters  of  6.5,  8.5  and  9.0  bohr.  It  is  seen  from 
Fig.  13  that  the  energy  profile  is  well  captured  using  the 
technique  within  expected  error  bounds  as  discussed  later 
in  this  section.  The  significant  advantage  of  using  inter¬ 
polated  energies  is  that  it  does  not  utilize  any  signifi¬ 
cant  computational  resource.  This  is  due  to  the  fact  that 
all  heavy  ab-initio  calculations  are  performed  beforehand 


FIG.  12:  The  three-body  cluster  configurational  space  for 
Platinum  colored  by  the  cohesive  energies  (EC(X i,  A'2,  A'3)  = 
E*( A'i,  A'2,  A'3)  — 53?=1  E*(Xi))  computed  from  ab-initio  sim¬ 
ulations. 


FIG.  13:  Comparison  of  variation  of  energies  with  lattice  pa¬ 
rameter  for  a  periodic  FCC  Pt  lattice  with  wMBE  calculations 
involving  cluster  energy  interpolation. 

and  the  data  is  stored  for  interpolation. 

f.  Analysis  of  accuracy  of  wMBE  with  interpolated  ab- 
initio  energies 

The  main  sources  of  error  in  the  wMBE  procedure 
are  the  errors  involved  in  the  interpolation  of  energies 
from  the  database,  fitting  weighting  coefficients  and  con¬ 
vergence  accuracy  for  the  ab-initio  energy  data.  The 
maximum  interpolation  error  over  any  element  (includ¬ 
ing  higher  dimensional  elements)  is  tightly  bounded  by 
Ctrmci  where  the  absolute  curvature  of  the  true  ab-initio 
energy  surface  is  bounded  in  each  element  f  by  a  con¬ 
stant  2 ct  and  rmc  is  the  minimum  containment  radius  of 
an  element.  For  the  2-atom  Pt  cluster  energy  surface, 
the  maximum  interpolation  error  was  0.03  mRyd.  Al¬ 
though  this  error  cannot  be  completely  eliminated,  we 
use  smaller  element  sizes  in  the  regions  where  large  en¬ 
ergy  variations  are  expected  in  order  to  reduce  the  inter¬ 
polation  error  for  larger  clusters.  The  convergence  accu¬ 
racy  of  self-consistent  DFT  calculations  of  small  clusters 
was  within  0.01  mRyd  in  all  cases.  In  order  to  study  the 
error  in  fitting  MBE  weights,  we  carried  out  a  leave-one- 
out  cross-validation  (CV)  procedure.  Here,  the  error  in 
reproduction  of  energies  is  studied  by  fitting  the  energy 
with  N  —  1  clusters  and  computing  the  error  in  reproduc¬ 
tion  of  energy  of  the  left-out  cluster  (£)).  The  process  is 


repeated  with  every  single  cluster  used  once  as  a  left-out 
cluster.  The  CV  error  is  computed  as  the  mean  error 
-k  \Etme  —  Ei\.  Compared  to  statistical  estimates 
such  as  variance,  CV  error  provides  a  more  reliable  esti¬ 
mate  of  future  performance  of  wMBE  when  energies  of 
new  clusters  need  to  be  predicted. 

Ab-initio  energies  of  N  =  300  randomly  generated  24- 
atom  Pt  clusters  were  used  for  testing  the  accuracy  of  the 
wMBE  procedure.  The  24  atoms  were  randomly  placed 
at  grid  points  spaced  7  bohr  apart  in  each  direction  over  a 
cube  of  105  bohr  length  and  ab-initio  energies  ( E )  for  use 
in  Eq.  (6)  are  computed.  The  mean  CV  error  for  third- 
order  MBE  was  found  to  be  0.381  Ryd  (15.9  mRyd  per 
atom).  The  mean  CV  error  during  cross  validation  for 
fourth-order  expansion  reduces  to  0.121  Ryd  (5.04  mRyd 
per  atom) .  The  average  cohesive  energy  per  atom  for  the 
complete  data  set  was  312.4  mRyd.  This  demonstrates 
convergence  towards  ab-initio  energies,  although  the  er¬ 
ror  may  still  be  significant  for  modeling  phenomena  such 
as  phase  transformations  where  accuracy  in  the  order  of 
mRyd  may  be  required. 

5.  Convergence  of  wMBE  for  a  binary  system  (a- 
alumina  AI2O3):  A  multi-body  expansion  is  constructed 
for  a-Alumina  ( AI2O3 )  system  using  cluster  energies 
computed  using  the  Streitz-Mintmire  (SM)  model25. 
Streitz-Mintmire  potential  is  a  many-body  functional 
that  merges  electrostatic  potential  with  an  embedded- 
atom  potential  to  describe  metal-oxide  energies,  a- 
Alumina  ( AI2O3 )  has  a  rhombohedral  primitive  unit  cell 
and  is  described  in  space  group  R3c  (no.  167)  with  two 
lattice  parameters  a,  b.  The  lattice  parameter  a  is  varied 
while  b  is  fixed  at  0.4856  bohr.  Figure  14  plots  the  varia¬ 
tion  of  energies,  computed  using  wMBE,  as  a  function  of 
lattice  parameter  a  for  a  2  x  2  x  1  cluster  of  a- Alumina. 
The  true  energies  as  computed  by  the  SM  model  at  each 
lattice  parameter  are  also  shown.  Four  energies  at  lat¬ 
tice  parameters  a  =  7.0,  7.2,  7.4  and  7.6  bohr  were  used 
to  compute  the  MBE  coefficients.  Within  this  range  of 
lattice  parameters,  energies  increase  linearly  as  indicated 
in  Fig.  14.  In  spite  of  this,  a  fourth-order  expansion  is 
able  to  represent  the  curvature  of  the  a- Alumina  energy 
profile  predicted  by  the  Streitz-Mintmire  (SM)  model.  In 
contrast  to  the  FCC  Pt  case  in  Fig.  7,  predicted  energies 
from  the  3rd  order  multi-body  expansion  is  not  able  to 
predict  the  energy  minima,  while  the  fourth-order  pre¬ 
dicts  the  lattice  parameter  a  =  6.6  bohr  accurately.  In¬ 
stead  of  the  SM  model,  ab-initio  calculations  of  isolated 
cluster  energies  ( E *)  could  have  been  used.  Since  we 
compute  energies  of  isolated  clusters  by  approximating 
a  periodic  lattice,  care  must  be  taken  to  avoid  the  in¬ 
fluence  of  lattice  Coulomb  potential  on  the  ionic  Al-0 
cluster  (due  to  finite  size  effects)  by  using  a  large  enough 
unit  cell.  DFT  calculations  were  avoided  in  this  example 
due  to  the  computational  complexity  of  handling  a  large 
number  of  plane  waves  because  of  the  sharply  peaked  va¬ 
lence  states  in  oxygen  and  requirement  of  a  large  unit 
cell.  wMBE  of  a  binary  metallic  system  that  uses  ab- 
initio  calculations  is  reported  in  the  next  example. 


FIG.  14:  Comparison  of  variation  of  energies  with  lattice  pa¬ 
rameter  for  a  2  x  2  x  1  supercell  of  a-alumina  (space  group 
R3c)  using  3rd  and  4th  order  wMBE.  The  true  energies  and 
the  energies  used  for  computing  MBE  coefficients  are  indi¬ 
cated. 


6.  wMBE  of  a  binary  (Au-Cu)  system  using  inter¬ 
polated  ab-initio  energies :  This  example  demonstrates 
structure  optimization  to  find  the  lattice  constants  for 
FCC  CuAii3  system  (space  group  Pm3m,  no.  221)  us¬ 
ing  interpolated  energies  of  clusters  computed  from  first 
principles  DFT  calculations.  As  in  the  case  of  FCC  Pt, 
the  energies  E  used  for  fitting  Eq.  (5)  are  obtained  from 
self-consistent  DFT  calculations  of  a  periodic  unit  cell.  A 
6x6x6  (864-atom)  FCC  cluster  is  considered  to  approx¬ 
imate  the  periodic  lattice,  and  MBE  expansion  is  con¬ 
structed  using  energies  interpolated  from  the  tessellated 
configurational  space.  As  an  example,  the  cohesive  en¬ 
ergy  ( Ec )  variation  with  cluster  specifiers  {R12,  R23,  R13) 
in  the  configurational  space  for  3-atom  Cu  —  Cu  —  Au 
and  Cu  —  Au  —  Au  clusters  is  shown  in  Fig.  15(a)  and 
(b),  respectively.  Apart  from  the  9  constraints  discussed 
in  section  B,  the  inequalities  R23  <  R13  and  R12  <  R13, 
respectively,  are  additionally  used  to  account  for  cluster 
symmetries  in  the  space  shown  in  Fig.  15(a)  and  (b). 
The  lower  and  upper  cutoffs  used  for  constructing  these 
spaces  were  2.19  bohr  and  10.95  bohr,  respectively.  For 
single  atom-type  clusters  of  copper  or  gold,  the  upper  and 
lower  cutoffs  were  fixed  at  0.3  and  1.5  times  the  lattice 
parameters  of  pure  FCC  Cu  and  Au  lattices. 

Figure  16  shows  comparison  of  the  energies  computed 
using  3rd  and  4th  order  wMBE  with  the  true  energies. 
Coefficients  in  the  multi-body  expansion  were  generated 
using  three  ab-initio  energy  calculations  of  a  periodic 
FCC  CuAus  lattice  with  lattice  parameters  of  8.6, 
8.7  and  8.8  bohr.  Similar  to  the  AI2O3  case,  the  3rd 
order  multi-body  expansion  is  not  able  to  capture  the 
energy  profile  of  FCC  CuAu3  whereas  a  fourth-order 
expansion  provides  a  reasonable  approximation  of  the 
energy  profile.  wMBE  approach  allows  computation  of 
the  energy  of  large  systems  with  accuracy  subject  to  the 
errors  discussed  previously.  Cross-validation  accuracy 
for  this  system  using  a  similar  procedure  as  described 
before  was  also  carried  out.  We  employed  300  random 
clusters  of  24-atom  CuAu3  clusters  for  testing  the 
accuracy  of  the  wMBE  procedure.  The  24  atoms  were 
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pansions  (MBE)  were  improved  by  weighting  ener¬ 
gies  obtained  from  various  orders  of  atom  interac¬ 
tions  in  a  new  method  called  weighted  multi-body 
expansion  (wMBE). 

•  In  contrast  to  methods  such  as  cluster  expan¬ 
sion  that  involve  the  ordering  degrees  of  freedom, 
wMBE  focuses  on  the  positional  degrees  of  free¬ 
dom.  This  allows  one  to  explicitly  model  structural 
relaxations. 


FIG.  15:  The  three-body  cluster  configurational  space  for  (a) 
Cu  —  Cu  —  Au  and  (b)  Cu  —  Au  —  Au  colored  by  the  cohe¬ 
sive  energies  (as  defined  in  Fig.  12)  computed  from  ab-initio 
simulations.  The  two  spaces  have  different  geometries  due  to 
different  underlying  symmetries. 


FIG.  16:  Comparison  of  variation  of  energies  with  lattice  pa¬ 
rameter  for  a  periodic  FCC  CuAus  lattice  with  wMBE  cal¬ 
culations  involving  cluster  energy  interpolation. 


•  Database  interpolation  techniques  are  demon¬ 
strated  for  accelerating  computation  of  energies  us¬ 
ing  multi-body  expansions.  For  the  first  time  in 
literature,  multi-body  expansions  were  computed 
directly  from  ab-initio  energies  of  small  clusters 
to  model  energies  of  Platinum  and  a  binary  al¬ 
loy  (Au-Cu)  system.  The  quality  of  the  expansion 
was  quantified  using  leave-one-out  cross-validation 
technique. 

•  The  technique  involves  considerably  lesser  compu¬ 
tational  cost,  with  no  requirement  of  periodicity, 
and  hence,  could  be  used  to  perform  more  accu¬ 
rate  NVE  or  NVT  molecular  simulations  of  metallic 
clusters  and  complex  phase  structures  compared  to 
other  commonly  used  position-dependent  potential 
approximations.  We  are  currently  working  on  data- 
adaptive  hierarchical  interpolation  which  would  al¬ 
low  us  to  build  higher  (5+)  order  potentials  that 
would  lead  to  improved  accuracy. 


randomly  placed  at  grid  points  spaced  7.5  bohr  apart  in 
each  direction  over  a  cube  of  112.5  bohr  length.  A  cross 
validation  error  of  0.187  Ryd  (7.8  rnRyd  per  atom)  was 
achieved  when  a  fourth-order  expansion  was  used.  The 
average  cohesive  energy  per  atom  for  the  complete  data 
set  was  207.1  mRyd.  This  error  may  be  significant  for 
modeling  phenomena  such  as  phase  transformations,  but 
wMBE  is  a  good  replacement  for  empirical  potentials 
in  several  other  multiscale  modeling  applications  where 
reasonable  accuracy  is  required.  Although  use  of  higher 
(5+)  body  interactions  is  expected  to  improve  the  fit,  it 
greatly  increases  computational  overhead  in  tessellation 
and  data  generation.  We  are  currently  working  on  the 
use  of  data-adaptive  hierarchical  interpolation  to  address 
this  issue. 


E.  Conclusions 

Developments  presented  here  advance  the  existing 
state  of  the  art  in  multi-body  expansion  technique  for 
representation  of  energies  of  alloy  systems  through  the 
following  new  contributions: 

•  The  convergence  characteristics  of  multi-body  ex- 


APPENDIX  A:  AB-INITIO  CALCULATIONS 

Ab  initio  electronic-structure  calculations  were  car¬ 
ried  out  using  density  functional  theory  in  the  local 
density  approximation,  as  implemented  in  the  PWscf 
package,  using  Perdew-Zunger  parameterization  of  the 
exchange  correlation  energy  and  Rabe-Rappe-Kaxiras- 
Joannopoulos26  (ultrasoft)  pseudopotential.  Kohn-Sham 
orbitals  were  expanded  in  a  plane  wave  basis  up  to  an 
energy  cutoff  calculated  to  ensure  convergence.  Brillouin 
zone  integrations  were  carried  out  using  single  fc-point 
calculation  and  Methfessel-Paxton  first-order  spread¬ 
ing27.  The  cell  size  is  taken  to  be  sufficiently  large  to 
effectively  simulate  an  isolated  cluster.  For  Platinum, 
the  energy  cutoff  was  244.8  eV,  the  cell  size  was  taken  as 
four  times  the  maximum  size  of  the  cluster. 
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Abstract 

A  Bayesian  scheme  to  fit  Potential  Energy  Surface  of  clusters  of  N  atoms  is 
proposed  using  a  permutationally  invariant  polynomial  basis.  The  evidence 
approximation  is  employed  to  fit  the  missing  prior  parameters  and  identify 
the  length  scale.  Distance  geometry  techniques  are  introduced  to  efficiently 
sample  the  configuration  space.  The  Bayesian  variance  is  used  to  quantify 
the  informational  content  of  each  point  in  the  configuration  space  leading  to 
an  efficient  adaptive  scheme  that  minimizes  the  required  number  of  expensive 
ab  initio  calculations.  Objective  stopping  criteria  are  provided. 

Keywords:  PES  interpolation,  model  selection,  invariant  polynomial  basis 


1.  Introduction 

The  potential  energy  surface  (PES)  plays  a  central  role  in  the  computa¬ 
tional  simulation  of  all  types  of  atomic  interactions  of  interest.  Once  the  PES 
is  constructed,  Molecular  Dynamics  (MD)  or  Monte  Carlo  (MC)  methods  can 
be  employed  to  investigate  the  system’s  dynamical  behavior.  Of  course,  the 
accuracy  of  such  simulations  depends  crucially  on  the  accuracy  of  the  PES 
used.  The  ideal  PES  is  the  so-called  Born-Oppenheimer  PES  obtained  by 
solving  the  Schrodinger  equation  using  the  adiabatic  approximation  [1],  Di- 
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rectly  using  the  ab  initio  PES  in  simulation  is  computationally  infeasible  for 
all  but  extremely  simple  systems. 

As  the  number  of  atoms  in  the  system  increases,  simple  functional  forms 
based  on  physical  models  are  a  necessary  choice.  Their  parameters  are  fitted 
using  a  small  set  of  experimental  data:  bond  energies,  bond  distances  and 
angles,  clastic  moduli,  vibrational  frequencies  etc.  Such  PES  give  qualitative 
descriptions  of  the  system  and  their  applicability  depends  closely  on  what 
range  of  experimental  data  was  used  to  fit  their  parameters,  i.e.  they  are  not 
transferable. 

Recently  it  has  become  possible  to  obtain  the  PES  directly  from  ab  ini¬ 
tio  calculations  for  relatively  small  and  not  too  complex  systems.  The  basic 
problems  addressed  in  the  literature  consist  of  1)  determining  the  important 
areas  of  the  configuration  space,  2)  Ending  the  minimum  number  of  the  con¬ 
figuration  space  points  required  to  obtain  an  accurate  PES  and  3)  fitting 
the  electronic  energies  to  an  analytic  model.  The  first  problem  is  usually  ad¬ 
dressed  by  employing  classical  trajectories  whose  initial  points  are  selected  to 
match  the  distribution  of  these  variables  under  the  conditions  of  the  experi¬ 
ments  being  investigated  [2],  Ab  initio  calculations  are  performed  on  a  set  of 
system  configurations  along  these  trajectories.  Alternatively,  an  approximate 
empirical  PES  can  be  used  to  initiate  the  trajectory  sampling  of  the  config¬ 
uration  space  [3].  The  second  problem,  testing  the  convergence  of  the  PES, 
is  performed  by  computing  various  dynamical  properties  of  the  system  and 
examining  their  invariancy  with  respect  to  the  database  size  [4].  Finally  sev¬ 
eral  methods  have  been  proposed  to  accurately  fit  ab  initio  databases.  Many 
methods  assume  parametrized  analytical  forms  for  the  surface.  Such  are  the 
many-body  expansion  method  [5]  and  the  recently  successful  multinomial  ex¬ 
pansion  method  [6,  7,  8].  Other,  basis  free,  methods  are  a)  moving  Shepard 
interpolation  techniques  [4,  9,  10,  11,  12,  13],  b)  reproducing  kernel  Hilbert 
spaces  [14],  c)  interpolating  moving  least  squares  [15,  16,  17,  18,  19,  20,  21] 
and  d)  Neural  Networks  methods  [22,  23,  24,  25]. 

Despite  the  appeal  and  usefulness  of  the  above  mentioned  techniques 
they  all  suffer  from  the  curse  of  dimensionality:  it  is  practically  impossible 
to  construct  a  PES  for  a  multi-atom  system.  To  overcome  this  barrier  we 
employ  the  Multi-Body  Expansion  (MBE)  technique  [26].  MBE  provides  a 
systematic  framework  in  which  the  total  energy  of  a  multi-atom  system  is 
represented  as  a  summation  over  potentials  of  isolated  clusters,  with  series 
terms  involving  pair,  three-body,  four-body,. . . ,  IV-body  potentials.  This  re¬ 
sults  in  structure  independent,  fully-transferable  many-body  potentials  [27]. 
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The  iV-body  potentials  can  be  constructed  using  the  Mobius  transformation 
from  the  A'-body  PES,  where  K  =  1,2, ... ,  N.  ffowever,  building  the  N- 
body  potentials  is  not  an  easy  problem:  1)  The  order  of  the  expansion  is  a 
priori  unknown  and  although  it  is  small  for  rare-gases,  it  is  of  relatively  high 
order  in  metals  [28].  2)  The  accuracy  of  the  IV-body  potential  depends  in  a 
complicated  manner  on  the  accuracy  of  the  A'-body  PES  for  A'  =  1, ...  ,N. 
Both  these  problems  pose  daunting  restrictions  on  the  applicability  of  the 
method  since  they  require  a  very  accurate  PES  fitting  scheme  that  utilizes 
as  little  as  possible  electronic  structure  calculations.  Furthermore,  for  the 
constructed  potentials  to  be  of  any  use,  the  analytical  form  of  the  PES  needs 
to  be  able  to  capture  a  wide  range  of  the  configuration  space.  Consequently, 
there  have  been  no  published  reports  of  a  multibody  expansion  constructed 
explicitly  from  first-principle  calculations.  The  main  goal  of  this  work  is  to 
address  exactly  these  issues. 

To  this  end,  we  propose  a  Bayesian  variant  of  the  multinomial  expansion 
method  for  the  construction  of  the  A'-body  PES.  The  multinomial  expan¬ 
sion  method  provides  an  analytical  functional  form  for  the  PES  satisfying 
permutation  invariance  with  respect  to  like  atoms.  Most  schemes  can  ac¬ 
count  for  permutation  invariance  only  by  explicitly  replicating  all  possible 
permutations  of  each  data  point.  As  a  result,  permutation  invariance  of 
like  atoms  is  learned  from  the  data  despite  the  fact  that  it  constitutes  a 
well-known  property  of  any  PES  -  a  well-established  prior  knowledge.  The 
accuracy  of  the  fitting  procedure  is  further  enhanced  by  introducing  a  Lin¬ 
ear  Bayesian  Regression  scheme  and  the  evidence  approximation  to  optimize 
the  scale  parameter  of  the  Morse-variables.  This  avoids  overfitting  and  in¬ 
creases  the  predictive  capabilities  of  the  PES.  An  additional  benefit  of  the 
Bayesian  framework  is  that  it  provides  us  with  a  way  to  quantify  the  in¬ 
formational  content  of  each  point  of  the  configuration  space.  Our  lack  of 
knowledge  about  the  energy  value  of  a  particular  point  of  the  configuration 
space  is  associated  with  the  variance  of  the  Bayesian  prediction.  This  in¬ 
formation  is  then  used  to  adaptively  select  new  data  points  to  be  included 
into  the  fitting  procedure  until  objective  stopping  criteria  are  met.  Distance 
Geometry  techniques  ([29,  30,  31])  are  introduced  to  facilitate  the  sampling 
of  the  configuration  space.  In  our  numerical  examples,  we  demonstrate  that 
this  improved  multinomial  expansion  greatly  increases  the  accuracy  of  the 
obtained  PES  and  reduces  the  number  of  required  electronic  calculations. 

In  Section  2.1  we  introduce  the  distance  matrix  as  our  variables  of  choice 
to  describe  the  configuration  space  and  we  discuss  why  it  is  necessary  to 
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map  them  to  the  Morse  variables.  In  Section  2.2  we  discuss  the  multino¬ 
mial  expansion  method  and  give  a  simple  illustration  on  how  the  Monomial 
Symmetrization  Approach  (MSA)  [8]  can  be  used  to  construct  a  basis  of 
permutationally  invariant  polynomials.  In  Section  2.3  we  mathematically 
define  the  configuration  space  as  the  set  of  distance  matrices  satisfying  cer¬ 
tain  bounds  and  introduce  Distance  Geometry  techniques  that  allow  us  to 
sample  from  it.  Section  2.4  describes  the  Linear  Bayesian  Regression  scheme 
and  Section  2.5  employs  the  evidence  approximation  in  order  to  select  all 
the  missing  parameters  of  the  model.  In  Section  2.6  we  propose  objective 
measures  to  test  the  goodness  of  fit,  and  in  Section  2.7  we  show  how  one 
can  use  the  Bayesian  variance  to  adaptively  select  new  configuration  points. 
Finally  in  Section  2.8  we  introduce  the  MBE  methodology  and  show  how  our 
framework  can  be  applied  to  the  construction  of  the  iV-body  potentials  from 
the  A'-body  PES. 

2.  Methodology 

Consider  a  cluster  of  N  atoms  not  necessarily  of  the  same  type.  Atom 
i  can  be  described  by  its  cartesian  coordinates  ry  e  R  and  its  type  at.  In 
order  to  simplify  the  notation  we  will  refer  only  to  the  cartesian  coordinates 
R  =  (ri, . . . ,  rN)  but  keeping  always  in  mind  that  these  are  associated  with 
particular  atomic  types. 

2.1.  Choice  of  variables 

We  choose  to  describe  the  configuration  space  in  terms  of  interatomic 
distances.  This  is  not  the  optimal  choice  since  we  use  ^N(N  —  1)  variables 
instead  of  the  3N  —  6  independent  degrees  of  freedom  1 .  However,  interatomic 
distances  provide  a  nice  way  to  sample  the  configuration  space  using  Distance 
Geometry  techniques  (see  Section  2.3)  and  guarantee  translation  and  rotation 
invariance. 

Let  dij  be  the  distance  between  atoms  i  and  j,  i.e. 

dij=\ri-rj\2l  (1) 

where  |  ■  |2  is  the  Euclidean  norm  of  R3.  The  symmetric  matrix 

D  =  (di:i),  (2) 
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is  called  the  distance  matrix  corresponding  to  the  Cartesian  coordinates  R  = 
(r i, . . . ,  7\v).  Since  D  has  a  zero  diagonal,  it  contains  exactly  \N(N  —  1) 
elements. 

Onr  plan  is  to  expand  the  energy  surface  in  a  polynomial  basis  (see  Section 
2.2).  Unfortunately  using  d^  as  the  variables  of  the  polynomials  will  be  an 
unphysical  choice  since  as  atoms  i  and  j  move  far  apart  ( dij  — >  oo),  the 
energy  would  diverge.  For  this  reason,  we  introduce  one  more  coordinate 
change,  the  so  called  Morse  variables  exp(— A  dij).  To  simplify  the  notation 
let  K  =  \N(N  —  1)  and  define  the  K- dimensional  vector  z  —  (z1} ,  zk)  by 

zx  =  e~Xdl2 ,  r2  =  e~Xdl3  ,...,zK  =  e~XdK-^K .  (3) 

The  newly  introduced  parameter  A  is  related  to  the  scale  of  the  problem.  We 
consider  A  as  an  unknown  parameter  to  be  inferred  from  the  data. 


2.2.  Choice  of  basis  functions 

Following  the  recent  advances  in  the  held,  we  introduce  a  polynomial  basis 
on  the  2  variables  that  is  invariant  with  respect  to  permutation  of  atoms  of 
the  same  type.  In  [7]  computational  algebra  techniques  are  used  to  find  a 
basis  on  the  space  of  invariant  polynomials  using  commercial  computational 
algebra  software  like  MAGMA  [32],  In  this  work,  we  follow  the  Monomial 
Symmetrization  Approach  (MSA)  of  [8]  to  construct  the  polynomials. 

To  introduce  MSA  consider  a  cluster  of  4  atoms  of  the  same  type.  In  this 
case  the  variable  describing  the  system,  2,  is  6-dimensional.  The  energy  is 
approximated  by  a  sum  of  monomials  up  to  degree  k 


E(z) 


k 


ri  a  b  c  d  e  f 


a-\-b-\-c-\-d-\-€.-\-  f — 0 


(4) 


where  a,  6,  c,  d,  e  and  /  are  all  non-negative  integers  and  Ca^Ctdiej  is  the 
coefficient  of  the  monomial  z^z^z^zfz^zg .  This  expression  is  clearly  not  in¬ 
variant  with  respect  to  permutations  of  the  atoms  for  arbitrary  choices  of 
coefficients.  However,  if  we  demand  to  have  permutation  invariance,  it  turns 
out  that  many  of  these  coefficients  must  be  the  same.  To  give  a  concrete 
example,  suppose  we  permute  the  atoms  as 


(ffi,  r2,  r3,  r4)  ->  (r4,  r2,  rh  r3). 
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The  corresponding  permutation  in  the  distance  matrix  elements  is 

(yi2,  2/13,  2/14,  2/23,  2/24,  2/34)  ~ >  (2/24,  2/l4,  2/34,  2/12,  2/23,  2/l3) 

and  so  the  2  variables  permute  as 

(^1,  z2,  z3,  Z4,  z5,  Z6)  ->  (z5,  Z3,  Z6,  ^1,  Z4,  Z2)- 

We  want  E(z)  to  be  invariant  under  such  a  permutation,  i.e. 

E(zi,  z2,  z3,  z4,  z5,  z6)  =  E(z 5,  z3,  z6,  z  1,  z4,  z2), 

which  can  happen  only  the  monomials  z^z^zfz^zg  and  z^z\z^zfz%z{  have 
the  same  coefficient,  i.e. 


Cayb,c,d,e,f  Cd,f,b,e,a,c- 

Listing  all  permutations  will  give  us  exactly  which  coefficients  should  be 
equal.  The  corresponding  monomials  will  sum  to  form  a  single  basis  func¬ 
tion  invariant  with  respect  to  the  permutation  group.  We  call  this  proce¬ 
dure  symmetrization  and  we  denote  the  symmetrized  sum  of  monomials  by 
S[Z]_Z2Z^z4zIzq]  and  their  common  coefficient  by  -Da,&,c,d,e,/-  Linder  this  nota¬ 
tion  the  energy  is  written  as 

k 

E{z)  =  J2  Da,ltcMS[zlzb2zlztzlzl},  (5) 

a+b-\-c-\-d+e+f=0 

where  the  summation  is  over  exponents  a,  b,  c,  d,  e  and  /  that  give  a  unique 
S[Z]_Z2Z^z4zIzq\.  To  avoid  unnecessary  complications  we  denote  these  basis 
functions  by  (j>i(z),i  —  1 ...  M  and  write  the  energy  as 

M 

E(z;w)  =  ^^Wi<f>i(z),  (6) 

i= 1 

where  Wi,  i  =  1, . . . ,  M  are  the  corresponding  coefficients.  These  polynomials 


6 


up  to  degree  3  are: 


Mz) 
Mz) 
Mz ) 
Mz) 
Mz ) 

Mz) 
Mz ) 
Mz) 

4>io(z) 


Mz) 


=  i, 

=  Z3Z4  +  Z2Z5  +  ZiZq, 

—  (+3  +  £4)  (+2  +  Z5)  +  Zi(z2  +  Z3  +  Z4  +  Z5)  +  (^2  +  Z3  +  Z4  +  ^5)^6; 

2  1  2  1  2  1  2  1  2  1  2 

—  %  +  +  ^3  +  Z4  +  Z5  +  %  5 

=  Z3Z4(Z5  +  +  +  (+3-2-4  +  ^2^5  +  (^2  +  Z3  +  Z4  +  Z^Zq) 

+  ^2  (^3  (-24  +  Z5)  +  Z§(Z4  +  ^6))j 

=  Z1Z2Z4  +  Z-4Z3Z5  +  ^2^3^6  +  -2-4 +5  z6j 
=  Z1Z2Z3  +  Z1Z4Z5  +  Z2Z4Z6  +  Z3Z5Z6, 

=  Z3Z4(Z3  +  Z4)  +  22^5(^2  +  Z5)  +  Z\Zq{z\  +  Zq), 

—  zl(z2  +  Z3  +  Z4  +  Z5)  +  +  Z4  +  (Z3  +  ^4)^5) 

+  Zi(z%  +  £3  +  Z4  +  Zg)  +  (^3  +  Z4  +  ^5)^6  +  (23  +  Z4 

+  ^5)^6  +  ^(^3  +  £4  +  2%)  +  ^2(^3  +  ^4  +  -2^), 

=  Z$  +  4  +  4  +  z4+z5+  4- 


The  whole  procedure  can  be  carried  out  using  Zhen  Xie’s  code  which  can 
be  found  at  [33].  One  may  notice  from  the  basis  functions  given  above,  that 
many  monomials  of  the  z%  variables  appear  again  and  again.  For  example 
z\  appears  in  05  and  in  0iO.  This  fact  is  more  pronounced  in  polynomials 
of  higher  degree.  Calculating  each  basis  function  from  scratch  would  re¬ 
sult  in  repetitive  calculations  of  the  same  monomials  and  hence  it  would  be 
highly  impractical.  One  of  the  extremely  useful  features  of  the  above  men¬ 
tioned  program  is  that  it  uses  the  algorithms  described  in  [8]  to  break  down 
the  computation  of  the  basis  functions  in  reusable  parts.  This  significantly 
reduces  the  computational  time  needed  for  their  calculation.  We  have  devel¬ 
oped  a  Python  script  that  uses  the  output  of  this  program  to  produce  C++ 
and  Matlab  code. 


2.3.  Sampling  the  configuration  space 

We  manage  to  sample  the  configuration  space  using  Distance  Geometry 
techniques  which  have  been  applied  successfully  to  nuclear  magnetic  reso¬ 
nance  (NMR)  structural  determination  problems  [34],  For  a  review  of  the 
topic  and  further  information  on  the  algorithms  described  here  see  [29]. 
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Algorithm  1  USAMPLE(_L,  U):  Samples  a  distance  matrix  from  the  con¬ 
figuration  space  C(L,  U ). 

Require:  L,  U  and  unif()  (a  random  number  generator). 

{Construct  a  matrix  D  satisfying  the  bounds  and  the  triangle  inequality} 
for  i  —  1  to  N  —  1  do 
for  j  —  i  +  1  to  N  do 

(Sample  an  element  i,j  uniformly  within  the  bounds.} 
dij  kj  +  unif()  *  (uij  -  kj) 

{Update  bounds} 

'fJ'ij  dij 
I’ij  dij 

{Enforce  triangle  inequality} 

FLOYD(L,  U) 

end  for 
end  for 

{Project  D  to  the  closest  distance  matrix  D} 

D  <-  MME(D) 

if  D  is  not  within  the  initial  bounds  then 
Go  to  the  beginning  of  the  algorithm. 

end  if 
return  D 


We  start  by  defining  the  configuration  space.  Not  every  possible  config¬ 
uration  is  of  physical  interest  to  us  -  e.g.  situations  where  two  atoms  are 

very  close  together  or  very  far  apart.  Furthermore,  we  may  wish  to  restrict 

our  attention  to  regions  of  the  configuration  space  of  particular  interest,  like 
average  bond  lengths  or  angles  measured  in  experiments.  A  very  natural 
way  to  impose  such  constraints  to  the  configuration  space  is  provided  via  the 
Lower  Cut-Off  Matrix: 

L  =  (hi)  (7) 

and  the  LTpper  Cut-Off  Matrix: 

U  =  ( ).  (8) 

We  define  the  configuration  space  to  be  the  collection  of  distance  matrices 
that  lie  between  these  two  bounds,  i.e. 

C{L,  U)  :=  {D  :  D  =  (dij)  is  a  distance  matrix  s.t.  Uj  <  dt]  <  ut] }  .  (9) 


Algorithm  2  FLOYD(X,  U):  Update  lower  and  upper  bounds  of  the  dis¬ 
tance  matrix  by  enforcing  the  triangle  inequality. 

Require:  L,  U 

{Loop  over  all  pairs  of  atoms  N  times.} 
for  k  —  1  to  N  do 
for  i  =  1  to  N  —  1  do 
for  j  =  %  +  1  to  N  do 
if  Uij  >  Uik  +  Ukj  then 

Uij  *  'U'ik  Ukj 

end  if 

if  lij  / ,k  u kj  then 

lij  ‘  lik  Ujkj 

end  if 

if  kj  <  ljk  -  Uki  then 

l/j  ljk  Uki 

end  if 

if  lij  >  then 

print  "Bad  Bounds” 
end  if 
end  for 
end  for 
end  for 
return  L,  U 


This  is  the  space  we  wish  to  sample.  Figure  2.3  shows  how  the  three  dimen¬ 
sional  configuration  space  of  a  three  atom  cluster  looks  like. 

Sampling  a  distance  matrix  within  specific  bounds  is  not  a  trivial  task. 
A  distance  matrix  corresponding  to  a  cluster  of  TV-atoms  needs  to  satisfy 
N  different  types  of  inequalities:  the  triangle  inequality,  and  corresponding 
inequalities  involving  four  distances,  five  distances  and  so  on. 

For  this  reason  we  use  a  type  of  proposal-rejection  sampling  technique. 
We  first  sample  a  matrix  satisfying  the  bounds  and  the  triangle  inequality 
and  we  project  it  to  the  nearest  distance  matrix  in  a  least  squares  sense.  If 
the  resulting  matrix  satisfies  the  bounds  we  stop,  otherwise  we  repeat  the 
procedure.  This  is  outlined  in  Algorithm  1.  To  sample  the  matrix  satisfying 
the  bounds  and  the  triangle  inequality  we  iterate  over  each  (i,j)  pair  of 
atoms.  We  sample  uniformly  a  distance  dij  between  and  u^.  We  set 
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Algorithm  3  MME(L,  U):  It  projects  D  to  the  closest  distance  matrix  D 
in  a  least  square  sense.  For  details  see  [citation  needed]. 

Require:  D  =  (d,j) 

Calculate  the  TV-dimensional  vector  dcm. 
for  i  =  1  to  N  do 

J  1  79  1  j2 


dc.m  ,i 

end  for 


J_  \  '  ;2  _l_  \  \  ' 

N  2^1  j= 1  aij  N2  2-*ij= 1  2^fc= 


d2 

j+ 1  UA 


If  D  is  a  true  distance  matrix,  ricm  is  tin 
center  of  mass. 

Calculate  the  N  x  N  matrix  VF  = 
for  %  =  1  to  N  do 
for  j  =  1  to  N  do 


is  the  distance  of  each  atom  from  the 


end  for 


3  ~  dii) 


end  for 

If  D  is  a  true  distance  matrix,  then  W  would  be  the  metric  matrix  of  the 
cluster  i.e.  W  =  (r,:  •  rj). 

Compute  the  three  greatest  eigenvalues  of  W,  Ai  >  A2  >  A3  and  the 
corresponding  eigenvectors  w1,w2,w3. 

Define  a  matrix  X  of  size  N  x  3. 

Let  X  =  (xi,  x-2,x^)  (Column  view), 
for  i  =  1  to  3  do 

Xi  <-  V^W; 

end  for 

Now  each  row  of  X  =  {r\V2  ■  ■  ■  ^w)T  (row  view)  contains  a  cartesian  rep¬ 
resentation  of  the  closest  metric  matrix  to  W. 

Calculate  the  distance  matrix  D  associated  with  rq, . . . ,  r^v- 

return  D 


the  bounds  corresponding  to  this  distance  equal  to  dij  (lij  =  Uij  =  dij )  and 
finally  refine  all  other  upper  and  lower  bounds  so  that  the  triangle  inequality 
is  not  violated  using  the  FLOYD  algorithm  (see  Algorithm  2).  Finally  the 
resulting  matrix  is  projected  to  the  nearest  distance  matrix  by  employing  the 
Metric  Matrix  Embedding  algorithm  (see  Algorithm  3).  Figure  2.3  shows  the 
histograms  of  the  interatomic  distances  obtained  using  USAMPLE.  Notice 
that  the  sampling  is  not  uniform,  but  all  regions  of  the  configuration  space 
are  covered. 
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2-4-  Fitting  the  data 

Suppose  we  have  a  total  of  S  data  points  (z^\  E^)f=1.  To  avoid  lengthy 
notation  we  will  refer  to  these  data  collectively  as 

2>=(z«  £«)f=1.  (10) 


Figure  1:  A  view  of  the  configuration  space  for  a  cluster  of  3  atoms. 
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Figure  2:  Histograms  of  the  interatomic  distances  obtained  using  the  USAMPLE  algo¬ 
rithm.  The  minimum  distance  is  set  to  1.2  A  and  the  maximum  distance  to  13  A.  It  is 
obvious  that  the  sampling  is  not  uniform  but  any  point  of  the  configuration  space  has  the 
probability  of  occurring. 
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In  this  section  we  will  outline  how  we  fit  T>  to  the  model  described  in  Section 
2.2.  The  usual  approach  (see  [6])  is  to  fix  the  scale  parameter  A  (usually 
to  2  —  3  Bohr)  and  minimize  the  least  squares  error.  This  procedure  is 
straightforward  to  implement  but  has  several  drawbacks: 

1.  it  can  lead  in  severe  overfitting  if  many  basis  functions  are  used. 

2.  it  does  not  provide  a  quantification  of  the  uncertainty  when  doing  pre¬ 
dictions. 

3.  it  does  not  give  a  systematic  way  to  select  the  scale  parameter  A  or  test 
alternative  coordinate  transformations  than  the  Morse  variables  (e.x. 
Uij  =  1/ di:j  as  in  [10]  ). 

We  propose  the  use  a  Bayesian  Linear  Regression  scheme  in  order  to  cope 
with  precisely  these  problems.  We  use  the  generalized  linear  model  given 
in  Eq.  (6)  with  the  additional  assumption  of  an  additive  noise 

E{z)  =  E(z]  w)  +  aZ,  (11) 

where  Z  ~  A/"(0, 1)  and  cr2  is  the  variance  of  the  data.  Of  course,  our  ab 
initio  calculations  are  deterministic  and  hence  have  no  noise.  However,  in 
case  the  basis  functions  are  not  adequate  to  describe  the  energy  surface  the 
model  will  fit  data  with  noise.  Not  being  able  to  interpolate  between  data 
points  of  the  surface,  it  will  assume  that  their  variability  is  due  to  a  big  o. 
That  is,  the  value  of  o 2  will  be  an  indicator  of  how  good  the  fit  really  is. 
Under  this  assumption  the  likelihood  of  the  real  energy  E(z)  is 

p(E\z,  w,  a2)  =  A f(E(z;  w ),  a2).  (12) 

Finally,  we  pose  a  Gaussian  prior  distribution  over  the  weights 

p(w)  =  jV(/io,  S0),  (13) 

where  Af(/i o,  So)  is  the  multivariate  normal  distribution  with  mean  /Uq  and 
covariance  matrix  S0.  It  turns  out  [35]  that  the  posterior  distribution  of  the 
weights  given  the  data  is  also  Gaussian 

p{w\V)=M{p,Yl),  (14) 

with  mean  and  inverse  covariance  matrix  given  by 

l1  —  S  (S0  1/x0  +  cr  2<hT  E )  , 

S"1  =  E^  +  a"2#^, 
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(15) 

(16) 


where  $  is  the  so  called  S  x  M  design  matrix: 


*<* is,)T  J 


with 

4>(z)  =  (0i(z), . . .  Am{z))t 
and  the  observed  energy  vector 

E  =  (£(1),...,£(5))T. 

Selection  of  the  Prior.  We  propose  starting  with  an  isotropic  Gaussian  prior 

Mo  =  0,  (17) 

So  =  a.1,  (18) 

where  a  is  an  additional  unknown  parameter  which  will  be  inferred  from 

the  data  and  I  is  the  M  x  M  unit  matrix.  Notice  that  the  prior  becomes 
uninformative  as  a  — >  0  allowing  for  a  very  flexible  choice  of  weights. 

Predictive  Distribution.  The  predictive  distribution  at  a  new  point  z  can  be 
calculated  by  integrating  out  the  weights  using  their  posterior  distribution. 
Not  surprisingly  this  is  a  Gaussian  also: 

p(E\z,  V,  a,  a2)  =  /  p(E\z,  w,  cr2)p(w\'D,  a)dw  =  A/^/i7  <t2(z)), 

(19) 

where  the  variance  <r2( z)  is  given  by 

a2  [z)  =  +  <f>(yZ)Ti(j)(z) .  (20) 

A  point  estimate  of  E(z)  can  be  given  by  the  mean  pT(f>(z).  Notice  that 
we  automatically  get  a  quantification  of  the  uncertainty  of  our  prediction 
through  cr2(z).  This  will  be  used  in  Section  2.7  in  the  selection  of  new  data 
points. 
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2.5.  Bayesian  Model  Selection 

In  the  previous  section  we  have  shown  how  the  model  fits  the  data  for 
a  fixed  choice  of  a  and  cr2.  The  careful  reader  would  have  noticed  that  the 
prediction  is  also  implicitly  dependent  on  the  choice  of  the  scale  parameter 
A,  i.e. 

z  =  Z\(D), 

where  D  is  the  distance  matrix.  In  a  fully  Bayesian  scheme  we  would  have 
to  define  prior  distributions  on  all  those  three  parameters  and  then  integrate 
over  them  to  get  the  fully  Bayesian  predictive  distribution.  This  is  not  an 
easy  task  to  perform.  In  this  section  we  will  motivate  the  so  called  evi¬ 
dence  approximation  [36]  which,  under  special  assumptions,  can  provide  an 
alternative  way  to  solve  the  problem. 

Suppose  we  have  introduced  some  prior  distribution  p(a,  a2,  A)  on  (a,  a2,  A). 
By  Bayes  Theorem  the  posterior  distribution  is 

p(a,  (3,  X\V)  =  p(V\a ,  a2,  A )p(a,  a2,  A),  (21) 

and  the  fully  Bayesian  predictive  distribution  is 


p{E\D,V) 


p(E\z\(D),  w,  a2)p{w\V,  a)p(a,  (3,  X\V)dwdadcr2dX. 


If  the  posterior  is  sharply  peaked  around  (a,cr2,A)  it  can  be  treated  as  a  5 
function  and  hence  the  above  integral  may  be  approximated  by 

Intuitively,  this  would  be  the  case  if  the  basis  functions  we  are  using  can 
describe  the  true  energy  surface  and  if  we  have  sufficient  data  at  our  disposal. 
In  this  case  predictions  can  be  made  using  Eq.  (19)  at  (a,  a2,  A). 

Now  the  problem  becomes  to  maximize  the  posterior  of  the  hyperparam¬ 
eters  (a,  cr2,  A)  given  by  Eq.  (21).  This  is  not  a  trivial  task  either.  However, 
if  the  prior  p(a,  cr2,  A)  is  relatively  flat  (which  would  be  the  case  if  we  have 
no  specific  preference  about  these  parameters  ),  then  this  maximum  will  be 
effectively  the  maximum  of  the  likelihood  function  (see  Eq.  (21)  again): 

C(a,  a2,  A)  =  p(V\a,  cr2,  A). 
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This  function  is  called  the  marginal  likelihood  or  the  evidence  function.  Our 
model  selection  problem  becomes: 

(a,  a2,  A)  =  arg  max  C(a,a2,  A).  (22) 

(a, o'2, A) 

In  principle  it  can  be  solved  by  conjugate  gradient  methods.  Unfortunately 
this  would  require  the  derivative  of  the  basis  functions  with  respect  to  A  which 
is  a  very  involved  calculation.  In  our  numerical  experiments,  we  exploit  the 
fact  that  for  fixed  A  there  exist  a  robust  algorithm  (see  Chapter  3.5  of  [35]) 
to  solve 

(a,  a2)  =  arg  max  C(a,  cr2,  A), 

(a,  o'2) 

in  order  to  pose  the  problem  as 

A  =  arg  max  I  max  C{a,  a2,  A)  )  .  (23) 

A  \(a,cr2)  '  J 

The  optimization  over  A  can  be  carried  out  using  Brent’s  method  (see  [37]). 


2. 6.  Error  Evaluation 

A  concrete  way  to  account  for  the  error  of  the  approximation  is  to  leave 
out  of  the  fitting  procedure  some  samples  and  compare  the  prediction  of 
their  energy  with  the  true  value.  The  measure  we  propose  to  use  is  the  mean 
square  error  of  the  energy  per  atom 


MSE(Aest) 


Attest 


where  we  have  left  out  of  the  fitting  procedure  Stest  data  points: 


V 


test 


Stest 
2=1  • 


(24) 


An  alternative  measure  we  will  also  use  is  the  maximum  absolute  error  of 
the  energy  per  atom 


MABSE(Dtest) 


max 

i<l<5test 


E, 


(0 


test 


E{z, 


(0  N 

’test  J 


Both  of  them  are  objective  measures  of  the  predictive  ability  of  the  fit  we 
have  achieved.  Ideally,  one  would  like  to  see  MABSE(Dtest)  becoming  ap¬ 
proximately  the  same  as  MSE(X>test)-  This  would  mean  that  there  are  no 
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outlying  energies  and  that  our  fitting  is  uniform.  In  Section  2.7  we  will  use 
MSE(Aest)  to  define  a  stopping  criterion  of  our  scheme.  It  will  provide  us 
with  a  definite  way  to  judge  if  we  need  more  data  points  or  we  have  reached 
the  maximum  predictive  ability  of  the  chosen  basis. 

2. 7.  Adaptive  selection  of  data  points 

Ab  initio  data  is  expensive  so  we  cannot  afford  to  waste  any.  Inspired 
by  the  seminal  paper  of  Sacks  [38],  we  propose  a  simple  experimental  design 


Algorithm  4  BFED:  Fit  the  PES  of  a  given  cluster. 

Require:  1)  Cluster  type  XnYm,  2)  bounds  of  the  configuration  space 
( L,U ),  3)  maximum  polynomial  degree  dmax  4)  initial  number  of  data 
points  Sinit,  5)  number  of  test  samples  Stest,  6)  number  of  MC  samples 
whose  variance  is  tested  in  every  cycle  S'mc;  V  number  of  samples  that  are 
added  to  the  fitting  procedure  after  its  cycle  S'add- 

Generate  using  USAMPLE(H,  U)  Stest  samples  and  calculate  their  ener¬ 
gies.  We  denote  them  with 


V 


test 


(E 


(i) 

test ) 


Stest 
i=  1 


Generate  using  USAMPLE(L,  U)  Sjnit  samples  and  calculate  their  ener¬ 
gies.  We  denote  them  with 


V 


(E{'\ 


gfS)  ^^init 


loop 

Fit  T>  as  described  in  Section  2.4  and  Section  2.5. 
if  MSE  latest  doesn’t  change  significantly  from  its  previous  value  then 
return 
else 

Generate  Smc  configuration  point  samples  using  USAMPLE(.L,  U). 
Calculate  their  variance  using  Eq.  (20). 

Find  the  S'add  of  the  S'mc  samples  with  the  maximum  variance. 
Calculate  their  energy  and  add  them  to  V. 

end  if 
end  loop 
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scheme.  The  goal  is  to  extract  information  from  existing  data  that  would 
help  us  select  new  points  of  the  configuration  points  that  will  improve  our 
fit  by  keeping  the  ab  initio  calculations  to  a  minimum. 

We  have  already  mentioned  that  the  Bayesian  scheme  described  in  Sec¬ 
tion  2.4,  provides  a  natural  way  to  quantify  our  uncertainty  at  any  point 
of  the  configuration  space  z  through  the  predictive  variance  a2  (z)  defined 
in  Eq.  (20).  This  variance  represents  our  lack  of  knowledge  about  the  PES 
and  not  the  real  error.  It  corresponds  to  the  real  error  (qualitatively)  only 
to  the  extent  that  the  selected  basis  can  actually  describe  the  true  energy 
surface.  We  postulate  that  the  inclusion  of  new  configuration  points  with 
high  Bayesian  variance  will  improve  the  fitting  achieved  as  measured  by  the 
MSE(Aest)  more  than  the  merely  adding  random  points.  This  is  valid  to 
the  extent  that  the  basis  functions  we  use  can  actually  describe  the  PES 
under  consideration  correctly.  Our  numerical  experiments  have  shown  that 
this  is  the  case  for  the  permutationally  invariant  polynomial  basis  alongside 
with  the  Morse  variables,  if  a  sufficient  number  of  basis  functions  are  used. 
We  have  found  that  the  number  of  require  energy  calculations  is  reduced 
significantly. 

Associating  the  value  of  cr(z)  with  the  informational  content  of  the  point 
z,  a  natural  strategy  is  to  add  to  the  scheme  the  z’s  that  maximize  it. 
We  use  a  plain  Monte  Carlo  procedure  to  identify  the  important  points  of 
the  configuration  space.  The  proposed  algorithm,  called  Brute  Force  Experi¬ 
mental  Design  (BFED),  is  extremely  simple  to  implement  since  all  it  requires 
is  the  ability  to  sample  the  configuration  space  and  calculate  the  Bayesian 
variance.  One  starts  with  some  random  configuration  points,  calculates  the 
corresponding  energy  and  fits  them  using  the  scheme  described  in  Section 
2.4  and  Section  2.5.  The  MSE  is  evaluated  on  some  random  energies  left 
out  of  the  fitting  procedure  as  described  in  Section  2.6.  Then  we  generate 
many  configuration  points  within  the  specified  bounds  using  the  US  AMPLE 
algorithm  and  calculate  the  predictive  variance  of  each  one  of  them  using  the 
fitted  surface.  This  is  computationally  negligible  compared  to  the  ab  initio 
calculations.  Finally  we  select  the  ones  with  the  maximum  Bayesian  vari¬ 
ance,  calculate  their  energies,  add  them  to  the  scheme  and  refit  2 .  If  the  MSE 


2The  refitting  procedure  doesn’t  have  to  start  from  scratch.  One  can  initialize  the 
a,  /?  and  A  parameters  at  the  previously  obtained  optimum  values  resulting  in  improve 
convergence  the  optimization  schemes.  However,  the  cost  of  the  fitting  procedure  is  so 
small  compared  to  the  ab  initio  calculations  that  making  it  faster  would  have  no  observable 
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of  the  new  fit  stops  changing,  we  stop  otherwise  we  repeat  the  procedure. 
This  is  summarized  in  Algorithm  4. 

2.8.  Multi-Body  Expansion 

One  can  easily  imagine  that  the  energy  of  a  cluster  of  N  atoms  with 
Cartesian  coordinates  R  =  (r*i, . . . ,  rjv)  €  R,:iAr  can  be  written  as  an  infinite 
series 

N 

E(ru...,rN)  =  ^2E(ri)  +  ^V[2)(ri,rj)+Y^  V(3\ri,rj,rk)  +  . . . .  (25) 

i=  1  i<j  i<j<k 

The  hW’s  are  functions  of  K  atoms  and  are  called  multibody  potentials. 
One  can  think  of  V^K\rn, . . . ,  rlK)  as  the  energy  added  to  the 

system  due  to  interactions  between  K  atoms.  It  seams  reasonable  that  after 
a  certain  index  P,  interactions  between  P+1  atoms  are  unimportant,  so  that 

V(p+1)  «  0.  (26) 

If  this  is  the  case,  then  we  call  this  expansion  a  P-order  expansion.  We 
denote  the  total  energy  of  an  TV-atom  system  using  a  P-order  expansion  by 
EP  =  EP(vi, . . . ,  rN).  We  write 

p 

EP{ru  ...,rN)  =  Y^  E{K){ri,  ...,rN),  (27) 

K=  1 

where 

E(K\ru...,rN)  :=  v{K)(rh,  ■  ■  • ,  nj-  (28) 

Finally,  the  V^’s  can  be  readily  found  by  using  the  Mobius  inversion  ap¬ 
proach  from  number  theory  [27,  39].  Mobius  inversion  has  been  used  pre¬ 
viously  for  the  extraction  of  potentials  from  energy  data  in  [40,  41].  We 
have 

I\ 

V{K\rl,...,rN)  =  Y^(~l)K~L  E(rii:...,riK),  (29) 

L= 1  ii<-”<tK 

where  with  E  we  denote  the  true  energy  function.  Once  the  multibody 
potentials  have  been  constructed,  the  P-order  energy  Ep  of  an  TV-atom 
system  can  be  calculated  by  using  Eq.  (27). 


benefit. 
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3.  Numerical  Examples 

3.1.  Pt  clusters  using  EAM  energy  and  GULP 

In  the  first  set  of  examples  we  wish  to  test  the  convergence  properties  of 
the  suggested  scheme.  Due  to  the  computational  burden  of  ab  initio  calcula¬ 
tions,  exhaustive  tests  can  only  be  performed  using  an  empirical  potential.  In 
what  follows  we  consider  solely  clusters  of  Pt  with  Embedded- Atom  Method 
(EAM)  [42]  potential  energy  as  implemented  in  GULP  [43].  The  minimum 
distance  between  Pt  atoms  is  taken  to  be  2  A  and  the  maximum  distance  13 
A.  All  energies  are  in  eV. 

First,  we  investigate  the  properties  of  the  adaptive  addition  scheme  of  Sec¬ 
tion  2.7  and  demonstrate  that  it  greatly  improves  the  quality  of  the  sampled 
energies.  Then,  we  test  the  convergence  of  fitting  scheme  with  respect  to  the 
polynomial  degree  of  the  basis  and  finally  we  compute  the  EAM  Pt7  PES  to 
assess  the  applicability  of  our  work  to  relatively  big  clusters. 

Adaptive  selection  of  data  points.  Our  first  goal  is  to  demonstrate  that  the 
adaptive  addition  of  data  points  (described  in  Section  2.7)  has  a  positive 
effect  on  the  accuracy  of  the  fit.  To  test  this  claim,  we  choose  a  particular 
Pt  cluster  and: 

1.  Generate  a  certain  number  Stes t  of  random  test  points  and  calculate 
their  energies.  These  are  left  out  of  the  fitting  procedure  and  are  used 
only  for  the  evaluation  of  the  error. 

2.  We  generate  <5),^  initial  random  points  and  calculate  their  energies. 

3.  LIsing  the  same  S^-dt  initial  random  points,  we  run  the  BFED  Algorithm 
for  different  adaptive  strategies  (varying  the  number  Smc  °f  samples 
among  which  we  select  the  points  that  should  be  added  to  the  dataset). 

4.  For  each  strategy,  we  trace  how  the  two  error  estimates  MSE  and  AMSE 
(described  in  Section  2.6)  vary  as  a  function  of  the  number  of  data 
points  used  to  fit  the  PES. 

We  perform  this  test  on  Pt3  and  Pt4  clusters.  In  both  cases  we  add 
Sadd  =  100  data  points  at  each  step  of  the  BFED  Algorithm,  set  the  basis  to 
polynomials  of  up  to  degree  10  and  we  test  the  same  six  different  strategies 
of  adding  new  data  points: 

1.  No  Selection.  Skdd  totally  random  points  are  added  at  each  cycle,  i.e. 

Smc  =  S^dd  =  100. 

2.  Choose  100  among  200,  i.e.  Smc  =  200. 
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3.  Choose  100  among  500,  i.e  Smc  —  500. 

4.  Choose  100  among  1000,  i.e  Smc  =  1000. 

5.  Choose  100  among  10000,  i.e  Smc  =  10000. 

6.  Choose  100  among  20000,  i.e  Smc  =  20000. 

For  Pt3,  the  number  of  test  points  is  set  to  Stest  =  1000  and  the  number  of 
initial  random  points  to  »S'init,  =  200.  For  Pt4,  the  number  of  test  points  is 
set  to  S'test  =  3000  and  the  number  of  initial  random  points  to  Sinit  =  500. 
In  Figure  3  and  Figure  4  we  plot  the  evolution  MSE  and  AMSE  respectively 
for  Pt3.  Figure  5  and  Figure  6  depicts  the  same  quantities  for  Pt4. 

These  figures  show  clearly  that  the  naive  random  selection  technique  is 
inferior  to  the  suggested  scheme.  One  can  notice  two  underlying  properties 
of  the  choice  of  Smc'- 

1.  The  error  drops  faster  as  a  function  of  the  number  of  data  points,  when 
Smc  is  initially  increased  but 

2.  there  is  a  natural  limit  to  this  effect.  In  both  tests,  setting  Sms  to  a 
value  greater  that  1000  does  not  refine  the  errors  any  more. 

The  former  property  provides  good  evidence  that  the  scheme  indeed  adds 
informationally  rich  data  points.  The  latter,  puts  a  barrier  on  this  improved 
performance. 

Finally,  in  Figure  7  we  plot  the  histograms  of  1000  sampled  Pt4  energies 
that  resulted  from  two  of  the  strategies  discussed  above:  the  random  addi¬ 
tion  strategy  Smc  —  100  shown  in  Figure  7(a)  and  the  “Choose  100  among 
10000”  strategy  Smc  —  10000  shown  in  Figure  7(b).  Notice  that  the  adaptive 
strategy  yields  a  considerably  broader  distribution  of  energies,  covering  low 
energy  (stable)  configuration  points  as  well  as  high  energy  (unstable)  ones. 
This  effect,  i.e.  the  more  uniform  sampling  of  the  configuration  space,  is 
the  reason  why  the  proposed  adaptive  scheme  actually  has  this  effect  to  the 
observed  errors.  As  the  Smc  initially  increases,  provides  better  sampling  of 
the  configuration  space  as  seen  by  the  distribution  of  the  energies.  However, 
after  a  critical  value  of  Smc  the  “optimum”  energy  distribution  has  been 
reached  and  increasing  it  further  has  no  effect. 

Convergence  with  respect  to  polynomial  degree.  An  important  question  that 
we  want  to  pose  is  to  what  extent  does  the  choice  of  the  polynomial  degree 
effects  the  accuracy  of  the  final  PES.  There  are  basically  two  reasons  why 
one  should  care  about  this: 
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Figure  3:  Pt3:  evolution  of  MSE(2?test)  as  the  number  of  samples  increases. 


1.  If  the  polynomial  degree  is  low,  the  basis  has  low  expressivity  and  it 
might  not  be  able  to  capture  the  real  PES  and 

2.  increasing  the  polynomial  degree  arbitrarily  might  result  in  overfitting 
which  would  yield  a  PES  with  limited  predictive  capabilities. 


Figure  4:  Pt3:  evolution  of  MABSE(2?test)  as  the  number  of  samples  increases. 
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Figure  5:  Pt4:  evolution  of  MSE(2?test)  as  the  number  of  samples  increases. 


To  resolve  the  first  problem,  the  natural  strategy  is  to  use  higher  order  poly¬ 
nomials.  In  the  case  study  that  follows  it  is  clearly  demonstrated  that  the 
second  problem  (overfitting)  is  missing  from  our  scheme.  This  is  due  to  its 
Bayesian  nature.  We  consider  a  Pt4  cluster  and  fit  its  PES  using  polynomials 


Figure  6:  Pt4:  evolution  of  MABSE(2?test)  as  the  number  of  samples  increases. 
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(a)  (b) 

Figure  7:  Pt4:  The  histograms  shows  1000  data  points  used  in  the  fitting  procedure.  7(a) 
shows  data  points  selected  completely  at  random  while  7(b)  shows  the  end  result  of  the 
proposed  scheme  when  selecting  100  out  of  10000  data  points. 


of  4th,  6th,  7th,  8th  and  10th  degree.  For  each  degree  we  use  exactly  the 
same  Stest  =  100  and  5)^  =  1000  test  and  initial  data  points,  respectively. 
The  adaptive  addition  strategy  is  “choose  100  among  10000”  (S'add  =  100 
and  S'mc  =  10000).  To  give  a  pictorial  representation  of  the  convergence,  we 
fix  atoms  2,  3  and  4  to 


r2  =  (1-1, 1-1,0), 
r3  =  (-1.1, 1.1,0), 
r4  =  (1.1, -1.1,0), 

and  we  plot  the  energy  as  a  function  of  the  position  of  the  1st  atom  allowing 
it  to  move  on  the  z  =  0  plane,  i.e. 

n  =  (x,y,  2.1). 

We  choose  this  particular  set  up  because  for  ry  =  (0,  0,  2.1),  we  obtain  close 
to  the  tetrahedral  stable  configuration  of  Pt4,  thus  the  test  is  performated 
in  a  region  of  the  configuration  space  of  a  high  variability  in  energy.  Figure 
8  depicts  the  corresponding  cut  of  the  PES  for  each  polynomial  degree.  The 
EAM  energy  is  also  shown  for  reference.  One  can  clearly  notice  that  after 
degree  7  the  picture  stabilizes.  Figure  9  shows  the  absolute  error  the  6,  7,  8 
and  10th  degree  PES.  This  is  evaluated  by  comparing  the  predicted  energy 
at  each  point  with  the  true  EAM  value. 
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Figure  8:  Showing  convergence  of  the  fitted  PES  for  a  Pt4  cluster  with  respect  to  the 
polynomial  degree  used. 


Figure  9:  Showing  the  absolute  error  of  the  fitted  PES  for  Pt4  for  different  polynomial 
degrees. 
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Demonstrating  the  scheme ’s  accuracy  for  big  clusters.  As  yet  another  case 
study,  we  test  the  accuracy  of  the  scheme  used  for  relatively  big  clusters.  We 
fit  the  EAM  PES  for  Pt7  using  5th  degree  polynomials,  Stest  =  100,  S^d  = 
100  and  S'mc  =  1000.  We  add  data  points  until  subsequent  values  MSE  do 
not  differ  more  than  10-3.  The  total  number  of  “electronic”  calculations  is 
2000.  In  Figure  10,  we  compare  the  fitted  energy  to  the  true  EAM  energy  as 
a  function  of  the  x-y  position  of  the  first  atom 

ri  =  (x,y,  0), 

while  keeping  the  rest  fixed  to 

r2  =  (2, 2, 2.1), 

7-3  =  (2, -2, 2.1), 

7-4  =  (-2, 2, 2.1), 

7-5  =  (-2, -2, 2.1), 

7-6  =  (0,0, 2.1), 

r7  =  (0,0, -2.1). 

Again,  this  is  a  region  of  high  variability  in  energy.  We  observe  that  the  PES 
remains  fairly  accurate  even  for  this  relatively  low  polynomial  degree. 
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Figure  10:  Pt7:  Plots  of  the  true  EAM  energy,  the  prediction  of  our  scheme  and  the 
absolute  error  in  energy. 
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3.2.  Pt  clusters  using  DFT  calculations 

In  what  follows,  ab  initio  calculations  were  performed  using  density  func¬ 
tional  theory  (DFT)  in  the  local  density  approximation  (LDA)  as  imple¬ 
mented  in  the  PWSF  software  package  [44],  using  the  Perdew-Zunger  parametriza- 
tion  of  the  exchange  correlation  energy  and  Rabe-Rappe-Kaxiras- Joannopoulos1  [45] 
(ultrasoft)  pseudopotential.  The  cell  size  was  taken  to  be  sufficiently  large 
(4  times  the  maximum  interatomic  distance,  a.k.a.  0.05  to  0.2  nano-meters  ) 
to  effectively  simulate  an  isolated  cluster.  The  energy  cut-off  was  22  Ry  (« 

300  eV).  One  fc-point  was  used  for  the  plane  wave  basis.  The  above  param¬ 
eters  where  chosen  so  that  the  accuracy  of  the  DFT  calculations  was  within 
O.OleV/atom.  All  distances  are  in  Bohr  and  energies  in  Ryd  unless  otherwise 
specified. 

We  computed  the  DFT  PES  for  Pt2,  Pt3,  Pt4,  Pt5  and  Pt6  clusters. 

The  data  collection  procedure  was  parallelized  using  MPI  [46]  so  that  each 
CPU  core  computed  a  single  data  point.  In  all  computations  we  used  Smc  = 

10000  and  Aa dd  =  256.  The  minimum  cut-off  distance  was  set  to  3.7  Bohr 
for  all  clusters.  The  maximum  cut-off  was  set  to  14  Bohr  for  Pt2  and  10 
Bohr  for  all  others  cases.  We  stopped  the  BFED  Algorithm  either  when 
MSE  per  atom  was  less  than  le  —  03  Ryd,  or  when  its  change  between  two 
consecutive  iterations  was  less  than  le  —  05  Ryd.  Table  1  shows  the  detailed 
parameters  of  each  run  along  with  the  computational  resources  that  were 
required.  The  vast  majority  of  the  computational  time  was  spent  in  ab  initio 
calculations.  For  example,  the  Pt5  PES  required  approximately  39  hours 
using  256  CPU  cores  from  which  only  6  minutes  was  spent  to  actually  fit 
the  data.  The  computational  effort  put  in  our  fitting  scheme  can  safely 
be  neglected  in  comparison  to  the  cost  of  the  DFT  calculations.  In  the 
Pt6  calculation  we  were  forced  to  use  a  smaller  polynomial  basis  because 
of  the  computational  burden  involved.  This  resulted  in  a  considerably  less 
accurate  PES.  In  several  cases  -  for  example  when  the  random  clusters  were  in 
highly  unstable  configurations  -  the  DFT  self-consistent  calculations  failed  to 
converge  within  100  iterations.  The  configurations  for  which  this  happened 
were  removed  from  the  data  set. 

After  assisting  the  accuracy  of  the  fitted  surfaces,  we  use  them  to  com¬ 
pute  stable  structures  of  small  Pt  clusters  and  compare  the  results  with  the 
literature.  We  continue  by  investigating  the  decomposition  of  the  PES  in 
A'-body  potentials  Vk ■  Finally,  we  put  these  potentials  together  to  explore 
the  predictive  capabilities  of  raw  MBE. 
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Accuracy  of  the  fitting  procedure.  Tabic  2  shows  the  errors  of  each  of  the  fitted 
PES  as  measured  using  the  test  points  left  out  of  the  scheme.  We  show  both 
the  Mean  Square  Error  as  well  the  Maximum  Absolute  Error  (Section  2.6). 
Notice  that  the  Maximum  Absolute  Error  is  of  the  order  of  le-03  Ryd  (~ 
0.01  eV  )  per  atom  which  is  exactly  the  accuracy  of  our  DFT  calculations. 
Figure  11  gives  a  pictorial  representation  of  the  goodness  of  fit.  For  a  given 
cluster  we  plot  the  true  DFT  energy  versus  the  energy  prediction  for  each 
test  data  point.  The  horizontal  axis  of  the  plots  corresponds  to  the  DFT 
calculated  energy  and  the  vertical  axis  to  the  prediction  of  the  energy  using 
the  fitted  PES.  The  straight  line  is  the  y  =  x  line.  The  deviation  of  each  red 
point  from  this  line  represents  the  error  of  the  corresponding  test  point.  It  is 
apparent  that  the  predictions  become  more  noisy  as  the  cluster  size  increases. 
However,  the  important  feature  that  one  should  notice  is  that  the  error  is 
uniformly  bounded  over  the  energy  range. 

Using  the  fitted  PES  to  predict  stable  structures.  As  an  elementary  applica¬ 
tion,  we  use  the  fitted  PES  to  predict  stable  Pt  structures  of  up  to  6  atoms. 
The  optimization  is  performed  using  the  Simulated  Annealing  technique  [47] 
as  implemented  in  the  GSL  library  [48].  The  quantities  we  report  are  the 
average  bond  length,  the  symmetry  of  the  structure  and  the  binding  energy. 
The  binding  energy  Eb  is  defined  by 

Eb  —  A'atom  —  A’cluster/A^,  (30) 

where  Eatom  —  E{  —  —52.157 Ryd  is  the  one  atom  energy  (see  next  paragraph 
for  details  on  how  we  get  calculate  this),  Aciuster  is  the  energy  of  the  stable 
cluster  and  N  the  number  of  atoms.  Table  3  summarizes  our  results.  For  easy 
comparison  with  the  results  in  the  literature  lengths  are  reported  in  A  and 
energies  in  eV.  Overall,  we  obtain  the  same  stable  structures  as  the  ones  found 


Cluster 

Degree 

No.  Basis 

Sinit 

*Stest 

<Sfinal 

CPU  Cores 

Time 

Pt2 

10 

11 

512 

128 

512 

128 

03:05:16 

Pt3 

10 

67 

1024 

256 

1024 

128 

09:31:53 

Pt4 

8 

195 

1024 

256 

1279 

256 

17:11:06 

Pt5 

8 

580 

1024 

256 

3071 

256 

39:30:57 

Pt6 

5 

86 

1024 

256 

1280 

256 

17:12:31 

Table  1:  Computational  details  of  the  fitting  procedure  of  Pt  clusters. 


Cluster 

MSE(Dtest) 

MSE(Dtest)  /  Atom 

MABSE(Aest) 

MABSE(T’test) / Atom 

Pt2 

2.2e-03 

l.le-03 

4.9e-04 

2.4e-04 

Pt3 

1.2e-02 

3.9e-03 

5.1e-03 

1.7e-03 

Pt4 

2.9e-02 

7.3e-03 

8.0e-03 

2.0e-03 

Pt5 

4.9e-02 

9.7e-03 

1.6e-02 

3.1e-03 

Pt6 

1.4e-01 

2.3e-02 

2.9e-02 

4.8e-03 

Table  2:  The  definition  of  the  errors 


Notation 

Pt3 

Pt4 

Pt5 

Pt6-1 

Pt6-2 

Structure 

-  c 

o 

y 

o 

o  c 

w 

O 

0 

°  0 

0 

W 

«  c 

Q  o  c 

Bond  Len.  (A) 

2.43 

2.54 

2.50-2.59 

2.50-2.7 

2.43-2.57 

Bind  En.  (eV/atom) 

3.25 

3.63 

3.94 

4.19 

4.23 

Symmetry 

c2v 

cs 

Dsh 

C2v 

c2v 

Table  3:  Stable  structures  of  Pt  clusters  predicted  using  the  fitted  PES. 

in  [49].  The  bond  lengths  predicted  are  slightly  smaller  than  the  ones  found 
in  the  literature  while  the  binding  energies  greater.  These  differences  were 
expected,  since  in  [49]  they  employ  the  generalized  gradient  approximation 
(GGA)  instead  of  the  LDA  we  use  in  the  present  work.  However,  the  behavior 
of  both  the  binding  energy  as  well  as  the  average  bond  length  as  a  function 
of  the  cluster  size  is  consistent  with  the  [49]  results.  A  comparison  is  shown 
in  Figure  12. 

Decomposition  of  the  PES  in  potentials.  The  values  of  the  potentials  Vk  can 
be  approximated  through  the  Mobius  transformation  Eq.  (29)  by  using  the 
fitted  PES  in  place  of  the  real  one.  We  use  the  values  thus  obtained  to  fit 
the  potentials  using  exactly  the  same  regression  scheme  as  the  one  used  to 
fit  the  PES.  The  polynomial  basis  of  the  -order  potential  Vk  is  chosen  to 
be  the  same  as  the  one  used  to  fit  to  PtK  PES.  The  one-atom  energy  E\{r) 
is  equal  to  half  the  limiting  value  of  the  fitted  energy  E2{r)  of  two  atoms  as 
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the  interatomic  distance  goes  to  infinity,  i.e. 


Ex 


-  lim  E2{r)  ss  —52.157  Ryd. 

2  r— >oo 


-314.6  -314.5  -314.4  -314.3  -314.2  -314.1  -314  -313.9  -313.8  -313.7 


(e)  Pt6 


Figure  11:  Testing  the  goodness  of  fit  using  full  DFT  calculations  of  random  clusters.  For 
Pt2  we  use  128  energies  and  256  for  Pt3,  Pt4,  Pt5  and  Pt6. 
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This  is  half  the  amount  of  energy  required  to  disassociate  the  two  atom 
system.  Figure  13  shows  the  original  )  and  the  derived  potential  V^r). 
Visualizing  the  higher  order  potentials  is  best  achieved  by  showing  a  two- 
dimensional  piece  of  the  configuration  space.  In  Figure  14  we  consider  a  Pt3 
cluster  and  plot  E3,  V2  and  V3  as  a  function  of  the  x-y  position  of  the  first 
atom 

ri  =  (x,y,3.7), 

while  keeping  the  rest  fixed  to 

r2  =  (-2.2, 0,0), 
r3  =  (2.2, 0,0). 

Figure  15  shows  £4,  V2,V3  and  V4  as  a  function  of  the  x  —  y  position  of  the 
first  atom 

ri  =  (x,y,  3.7), 

while  keeping  the  rest  fixed  to 

r2  =  (-3.7, 0,0), 
r3  =  (3.7, 0,0), 
r4  =  (0, 3.2,0). 

Figure  16  shows  E5,  V2,V3,  V4  and  V5  as  a  function  of  the  x  —  y  position  of 
the  first  atom 

n  =  (x,y,  3.7), 


(a)  Binding  Energy  (ev/atom)  (b)  Bond  Length  (ev/atom) 

Figure  12:  Comparison  of  binding  energy  and  average  bond  lengths  obtained  using  our 
fitted  PES  for  Pt2,  Pt3,  Pt4,  Pt5  and  Pt6  with  the  results  found  in  the  literature  [49]. 
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while  keeping  the  rest  fixed  to 


r2  =  (-3.7, 0,0), 
r3  =  (3.7, 0,0), 
r4  =  (0, 3.2,0), 
r5  =  (0, -3.2,0). 

Finally,  Figure  16  shows  E5,  V2,  V3,  V4,  V5  and  Ve  as  a  function  of  the  x  —  y 
position  of  the  first  atom 

ri  =  {x,y,  3.7), 

while  keeping  the  rest  fixed  to 


r2  = 

(0,0,0), 

rs  = 

(-6.4,  0,0), 

r4  = 

(6.4,  0,0), 

r5  = 

(3.2,  6.4,0), 

r&  = 

(-3.2,  6.4,0) 

The  key  observations  one  can  make  out  of  these  figures  are: 

1.  The  sign  of  the  interatomic  potential  is  systematically  reversed  from 
one  order  to  the  other,  i.e.  V2  and  V4  are  negative,  while  V3  and  V5  are 
positive. 

2.  VK  is  indeed  becoming  less  and  less  important  as  K  increases. 

3.  The  importance  of  each  Vk  becomes  greater  in  regions  of  the  configu¬ 
ration  space  of  low  energy. 

MBE  approximation  to  ab  initio  energies.  In  what  follows,  we  denote  with 
MBE-K  the  Multi-Body  Expansion  of  order  K  given  by  Eq.  (27)  and  Eq. 
(28)  of  Section  2.8.  We  wish  to  investigate  to  what  extent  MBE-K  can 
accurately  predict  electronic  energies  of  Pt  clusters  with  number  of  atoms 
N  >  K .  Figure  17  shows  the  predictions  of  MBE-2,  MBE-3  and  MBE-4  on 
256  Pt4  test  points.  Figure  18  shows  how  the  surface  of  Figure  15  can  be 
approximated  by  MBE  of  order  2,  3  and  4.  Figure  19  shows  the  predictions 
of  MBE-3,  MBE-4  and  MBE-5  on  256  Pt5  test  points.  Figure  20  shows 
how  the  surface  of  Figure  16  can  be  approximated  by  MBE  of  order  3,  4 
and  5.  Figure  21  shows  the  predictions  of  MBE-4,  MBE-5  and  MBE-6  on 
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256  Pt5  test  points.  Figure  22  shows  how  the  surface  of  Figure  ??  can  be 
approximated  by  MBE  of  order  4,  5  and  6.  Based  on  these  results  one  can 
conclude  that: 

1.  MBE-3  can  give  a  qualitatively  consistent  picture  of  the  Pt4  PES  but 
for  ab  initio  accuracy  MBE-4  is  needed  (Figure  18). 

2.  MBE-3  completely  fails  to  give  a  good  Pt5  PES.  MBE-4  gives  a  qual¬ 
itatively  consistent  PES  but  MBE-5  is  needed  to  achieve  ab  initio  ac¬ 
curacy  (Figure  20). 

3.  MBE-4  captures  many  important  features  of  the  Pt6  PES,  MBE-5  per¬ 
forms  much  better  but  without  ab  initio  accuracy  (Figure  22). 

4.  In  Figure  22  we  notice  that  MBE-6  does  not  reproduce  the  exact  Pt6 
PES. 

It  is  evident  that  MBE-4  or  higher  is  required  to  get  a  qualitatively  cor¬ 
rect  PES  for  Pt  clusters.  However,  the  applicability  of  relatively  low  order 
MBE  to  bigger  clusters  requires  further  investigation.  The  final  point,  is  a 
demonstration  of  another  difficulty  that  needs  to  be  addressed.  The  reason 
why  MBE-6  fails  to  reproduce  the  exact  Pt6  PES  is  the  small  errors  intro¬ 
duced  to  the  potentials  Vk  during  their  fitting  procedure.  This  error  is  of  the 
same  order  as  the  accuracy  of  the  ab  initio  calculation  but  it  propagates  in 
a  complicated  manner  through  the  Multi-Body  Expansion  Eq.  (27)  possibly 
amplifying  itself.  Even  if  the  error  of  each  Vk  is  reduced,  application  of  raw 
MBE  to  sufficiently  large  clusters  could  potentially  amplify  it.  The  nature 
of  this  propagation  is  largely  unknown  and  further  research  is  required  to 
address  it  properly. 


(a)  Pt2  PES  (b)  V2 

Figure  13:  Pt2:  Extracting  the  V2  potential  from  the  PES. 
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Figure  14:  Pt3:  Decomposition  of  the  energy  of  a  Pt3  cluster. 

Using  MBE  to  predict  the  energies  of  larger  clusters.  In  this  very  last  exam¬ 
ple,  we  sample  some  Pt  clusters  with  number  of  atoms  N  =  7,  8  and  10,  cal¬ 
culate  their  ab  initio  energies  and  compare  them  to  MBE-K  K  =  2,  3, 4,  5,  6. 
Figure  23  plots  jointly  all  the  results.  Notice  that  the  main  contribution  to 
the  expansion  comes  from  the  one  energy  terms  and  that  potentials  add  cor¬ 
rections  of  the  order  of  1  —  2  Ryd.  The  accuracy  after  MBE-3  for  Pt7  and  Pt8 
of  the  order  of  0.05  Ryd  (0.6  eV)  per  atom  which  is  below  our  goal  of  0.01  eV 
per  atom.  Furthermore,  notice  that  despite  the  fact  that  consecutive  MBE 
approximations  fluctuate  about  the  correct  value  of  the  energy,  they  don’t 
seem  to  converge  (up  to  order  6).  As  a  matter  of  fact,  the  oscillations  are 
more  pronounced  for  the  largest  PtlO  cluster.  We  believe  that  this  is  a  clear 
demonstration  of  the  unaccounted  propagation  of  the  error  in  the  potentials 
we  first  observed  in  the  results  of  the  previous  paragraph. 
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Figure  15:  Pt4:  Decomposition  of  the  energy  of  a  Pt4  cluster. 


Figure  16:  Pt5:  Decomposition  of  the  energy  of  a  Pt5  cluster. 
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Figure  17:  Pt4:  Comparing  the  prediction  of  successive  MBE  approximations  on  256  test 
points. 
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Figure  18:  Pt4:  Prediction  of  successive  MBE  approximations. 
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Figure  19:  Pt5:  Comparing  the  prediction  of  successive  MBE  approximations  on  256  test 
points. 
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Figure  20:  Pt5:  Prediction  of  successive  MBE  approximations 
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Figure  21:  Pt6:  Comparing  the  prediction  of  successive  MBE  approximations  on  256  test 
points. 
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Figure  22:  Pt6:  Prediction  of  successive  MBE  approximations 


Figure  23:  MBE  predictions  of  3  random  clusters  Pt7,  Pt8  and  PtlO. 
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4.  Conclusions 


The  construction  of  ab  initio  accurate  Potential  Energy  Surfaces  consti¬ 
tutes  a  problem  of  extreme  complexity,  albeit  of  vital  importance  towards 
the  search  for  materials  of  extremal  properties.  The  computational  cost  of  ab 
initio  calculations  makes  mandatory  the  search  for  sophisticated  techniques 
to  tabulate  ab  initio  data.  The  Multi-Body  Expansion  provides  a  mathe¬ 
matically  rigorous  framework  that  allows  one  to  break  down  the  PES  of  an 
arbitrarily  sized  cluster  to  interatomic  potentials  of  relatively  small  order. 
These  potentials  are  connected  to  Potential  Energy  Surfaces  of  small  order 
through  the  Mobius  transformation.  In  this  work,  we  presented  a  set  of  com¬ 
putational  techniques  to  efficiently  solve  the  PES  interpolation  problem.  The 
potentials  can  be  interpolated  using  the  exact  same  scheme. 

We  defined  mathematically  the  configuration  space  and  provided  algo¬ 
rithms  derived  from  the  Distance  Geometry  literature  to  efficiently  sample 
it.  The  recently  developed  multinomial  expansion  method  was  called  to  pro¬ 
vide  a  polynomial  basis  invariant  with  respect  to  permutations  of  like-atoms. 
This  provided  us  with  candidate  functions  that  satisfied  all  invariance  prin¬ 
ciples  of  an  energy  surface.  We  introduced  a  Bayesian  Linear  Regression 
scheme  to  fit  the  weights  of  the  polynomial  basis  as  well  as  the  evidence 
approximation  to  optimize  the  scale  parameter.  The  variance  of  the  predic¬ 
tion  was  used  to  devise  a  simple,  yet  effective,  way  to  adaptively  add  data 
points.  We  demonstrated  that  this  new  technique  considerably  improves  the 
quality  of  the  samples  obtained  and  outperforms  the  random  selection  of 
data  points.  Furthermore,  it  minimizes  the  number  of  ab  init  calculations 
required  enabling  us  to  construct  the  ab  init  PES  of  Pt  clusters  of  up  to  6 
atoms  using  a  reasonable  amount  of  computational  resources.  The  ab  initio 
PES  was  used  to  find  the  stable  structures  of  small  Pt  clusters  with  the  aid  of 
Simulated  Annealing.  The  results  were  found  to  be  in  good  agreement  with 
the  literature.  The  constructed  Pt  PES  was  also  used  to  fit  the  interatomic 
potentials  up  to  order  6.  It  was  shown  that  those  become  less  and  less  impor¬ 
tant  as  their  order  increases,  albeit  slowly  in  low  energy  regions.  We  used  the 
potentials  to  investigate  the  performance  Multi-Body  Expansions  of  various 
order  for  Pt  clusters  of  up  to  10  atoms.  It  was  demonstrated  that  interac¬ 
tions  of  at  least  5  atoms  are  required  to  qualitatively  describe  Pt  clusters. 
Finally,  we  observed  that  the  error  introduced  during  the  fitting  procedure 
of  the  interatomic  potentials  propagates  in  a  complicated  manner  through 
the  Multi-Body  Expansion  formula  making  its  naive  application  to  big  clus- 
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ters  questionable.  It  is  the  object  of  our  current  research  to  investigate  the 
propagation  of  this  error  through  the  MBE  formula  and  design  effective  tech¬ 
niques  to  filter  it  out.  We  believe  that  such  filtering  schemes  have  to  be  case 
specific  (different  for  each  material)  and  should  utilize  further  physical  infor¬ 
mation.  This  problem  constitutes  the  final  obstacle  towards  the  construction 
of  fully  transferable  potential  energy  surfaces  using  the  MBE  framework. 
The  impact  of  such  reduced-order  fully  transferable  PES  is  expected  to  be 
significant  in  the  search  for  new  materials,  exploring  materials  and  materials 
surface  design,  optimizing  mechanical  and  thermophysical  properties,  etc. 
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